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Summary 

The  objective  of  this  research  is  to  develop  an  efficient  and  accurate  methodology  to  resolve 
flow  non-linearity  of  fluid- structural  interaction.  To  achieve  this  purpose,  a  numerical  strat¬ 
egy  to  apply  the  detached-eddy  simulation  (DES)  with  a  fully  coupled  fluid-structural  in¬ 
teraction  model  is  established  for  the  first  time.  The  following  novel  numerical  algorithms 
are  also  created:  a  general  sub-domain  boundary  mapping  procedure  for  parallel  compu¬ 
tation  to  reduce  wall  clock  simulation  time,  an  efficient  and  low  diffusion  E-CUSP  (LDE) 
scheme  used  as  a  Riemann  solver  to  resolve  discontinuities  with  minimal  numerical  dissi¬ 
pation,  and  an  implicit  high  order  accuracy  weighted  essentially  non-oscillatory  (WENO) 
scheme  to  capture  shock  waves. 

The  Detached-Eddy  Simulation  is  based  on  the  model  proposed  by  Spalart  in  1997. 
Near  solid  walls  within  wall  boundary  layers,  the  Reynolds  averaged  Navier-Stokes  (RANS) 
equations  are  solved.  Outside  of  the  wall  boundary  layers,  the  3D  filtered  compressible 
Navier-Stokes  equations  are  solved  based  on  large  eddy  simulation(LES).  The  Spalart- 
Allmaras  one  equation  turbulence  model  is  solved  to  provide  the  Reynolds  stresses  in  the 
RANS  region  and  the  subgrid  scale  stresses  in  the  LES  region. 

An  improved  5th  order  finite  differencing  weighted  essentially  non-oscillatory  (WENO) 
scheme  with  an  optimized  e  value  is  employed  for  the  inviscid  fluxes.  The  new  LDE 
scheme  used  with  the  WENO  scheme  is  able  to  capture  crisp  shock  profiles  and  exact 
contact  surfaces.  A  set  of  fully  conservative  4th  order  finite  central  differencing  schemes 
are  used  for  the  viscous  terms. 

The  3D  Navier-Stokes  equations  are  discretized  based  on  a  conservative  finite  differ¬ 
encing  scheme,  which  is  implemented  by  shifting  the  solution  points  half  grid  interval  in 
each  direction  on  the  computational  domain.  The  solution  points  are  hence  located  in  the 
center  of  the  grid  cells  in  the  computational  domain  (not  physical  domain).  This  makes 
it  possible  to  use  the  same  code  structure  as  a  2nd  order  finite  volume  method.  A  finite 
differencing  high  order  WENO  scheme  is  used  since  a  finite  differencing  WENO  scheme 
is  much  more  efficient  than  a  finite  volume  WENO  scheme. 

The  unfactored  line  Gauss-Seidel  relaxation  iteration  is  employed  for  time  marching. 
For  the  time  accurate  unsteady  simulation,  the  temporal  terms  are  discretized  using  the 
2nd  order  accuracy  backward  differencing.  A  pseudo  temporal  term  is  introduced  for  the 
unsteady  calculation  following  Jameson’s  method.  Within  each  physical  time  step,  the 
solution  is  iterated  until  converged  based  on  pseudo  time  step. 

A  general  sub-domain  boundary  mapping  procedure  is  developed  for  arbitrary  topol¬ 
ogy  multi -block  structured  grids  with  grid  points  matched  on  sub-domain  boundaries.  The 
interface  of  two  adjacent  blocks  is  uniquely  defined  according  to  each  local  mesh  index 
system  (MIS)  which  is  specified  independently.  A  pack/unpack  procedure  based  on  the 
definition  of  the  interface  is  developed  to  exchange  the  data  in  a  ID  array  to  minimize  data 
communication.  A  secure  send/receive  procedure  is  employed  to  remove  the  possibility  of 
blocked  communication  and  achieve  optimum  parallel  computation  efficiency.  Two  terms, 
“Order’  and  “Orientation” ,  are  introduced  as  the  logics  defining  the  relationship  of  adja¬ 
cent  blocks.  The  domain  partitioning  treatment  of  the  implicit  matrices  is  to  simply  discard 
the  comer  matrices  so  that  the  implicit  Gauss-Seidel  iteration  can  be  implemented  within 
each  subdomain.  This  general  sub-domain  boundary  mapping  procedure  is  demonstrated 


to  have  high  scalability. 

Extensive  numerical  experiments  are  conducted  to  test  the  performance  of  the  numer¬ 
ical  algorithms.  The  LDE  scheme  is  compared  with  the  Roe  scheme  for  their  behavior 
with  RANS  simulation.  Both  the  LDE  and  the  Roe  scheme  can  use  high  CFL  numbers  and 
achieve  high  convergence  rates  for  the  algebraic  Baldwin-Lomax  turbulence  model.  For  the 
Spalart-Allmaras  one  equation  turbulence  model,  the  extra  equation  changes  the  Jacobian 
of  the  Roe  scheme  and  weakens  the  diagonal  dominance.  It  reduces  the  maximum  CFL 
number  permitted  by  the  Roe  scheme  and  hence  decreases  the  convergence  rate.  The  LDE 
scheme  is  only  slightly  affected  by  the  extra  equation  and  maintains  high  CFL  number  and 
convergence  rate.  The  high  stability  and  convergence  rate  using  the  Spalart-Allmaras  one 
equation  turbulence  model  is  important  since  the  DES  uses  the  same  transport  equation  for 
the  turbulence  stresses  closure. 

The  RANS  simulation  with  the  Spalart-Allmaras  one  equation  turbulence  model  is  the 
foundation  for  DES  and  is  hence  validated  with  other  transonic  flows  including  a  2D  sub¬ 
sonic  flat  plate  turbulent  boundary  layer,  2D  transonic  inlet-diffuser,  2D  RAE2822  airfoil, 
3D  ONERA  M6  wing,  and  a  3D  transonic  duct  with  shock  boundary  layer  interaction.  The 
predicted  results  agree  very  well  with  the  experiments.  The  RANS  code  is  then  further 
used  to  study  the  slot  size  effect  of  a  co-flow  jet  (CFJ)  airfoil. 

The  DES  solver  with  fully  coupled  fluid- structural  interaction  methodology  is  validated 
with  vortex  induced  vibration  of  a  cylinder  and  a  transonic  forced  pitching  airfoil.  For  the 
cylinder,  the  laminar  Navier-Stokes  equations  are  solved  due  to  the  low  Reynolds  number. 
The  3D  effects  are  observed  in  both  stationary  and  oscillating  cylinder  simulation  because 
of  the  flow  separations  behind  the  cylinder.  For  the  transonic  forced  pitching  airfoil  DES 
computation,  there  is  no  flow  separation  in  the  flow  field.  The  DES  results  agree  well  with 
the  RANS  results.  These  two  cases  indicate  that  the  DES  is  more  effective  on  predicting 
flow  separation. 

The  DES  code  is  used  to  simulate  the  limited  cycle  oscillation  of  NLR7301  airfoil.  For 
the  cases  computed  in  this  research,  the  predicted  LCO  frequency,  amplitudes,  averaged 
lift  and  moment,  all  agree  excellently  with  the  experiment.  The  solutions  appear  to  have 
bifurcation  and  are  dependent  on  the  initial  perturbation.  The  developed  methodology  is 
able  to  capture  the  LCO  with  very  small  amplitudes  measured  in  the  experiment.  This  is 
attributed  to  the  high  order  low  diffusion  schemes,  fully  coupled  FSI  model,  and  the  tur¬ 
bulence  model  used.  This  research  appears  to  be  the  first  time  that  a  numerical  simulation 
of  LCO  matches  the  experiment.  The  DES  code  is  also  used  to  simulate  the  CFJ  airfoil  jet 
mixing  at  high  angle  of  attack. 

In  conclusion,  the  numerical  strategy  of  the  high  order  DES  with  fully  coupled  FSI 
model  and  parallel  computing  developed  in  this  research  is  demonstrated  to  have  high  ac¬ 
curacy,  robustness,  and  efficiency.  Future  work  to  further  maturate  the  methodology  is 
suggested. 
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Chapter  1 


Introduction 

1.1  Background 

Flow  induced  structural  vibration  is  one  of  the  most  challenging  problems  affecting  the 
military  and  commercial  aircraft.  Due  to  the  very  complicated  non-linear  flow-structure 
interaction  phenomena,  there  is  a  lack  of  high  fidelity  computational  tools  to  study  the 
basic  physics  and  to  accurately  predict  the  structural  failure  limits.  For  airframe,  there  are 
problems  such  as  transonic  flutter,  limit  cycle  oscillation,  buzz,  buffet,  etc.  For  propulsion 
turbomachinery,  there  are  problems  such  as  high  cycle  fatigue  caused  by  forced  response 
or  stall  flutter,  etc.  Helicopter  rotor  blades  constantly  work  under  the  vibration  induced  by 
blade  wake  and  tip  vortexes.  The  development  of  advanced  methodologies  to  accurately 
simulate  fluid-structural  interactions  will  have  broad  impact  on  improving  the  performance 
of  various  aircraft. 

The  difficult  issue  that  the  fluid-structural  interaction  (FSI)  community  currently  faces 
is  the  non-linearity  caused  by  both  fluid  and  structure  [1].  The  aerodynamic  non-linearity 
poses  more  challenge  than  the  structural  one  [2,3].  These  problems  are  often  accompanied 
or  caused  by  the  complicated  flow  phenomena  such  as  shock  wave/turbulent  boundary  layer 
interaction,  flow  separation,  vortex  shedding,  etc. 

The  first  challenge  to  resolve  aerodynamic  non-linearity  such  as  shock  wave/turbulent 
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boundary  layer  interaction  and  flow  separation  is  the  turbulence  simulation.  The  second 
issue  is  the  high  accuracy  requirement  of  the  numerical  simulation  with  low  numerical  dif¬ 
fusion.  A  large  numerical  diffusion  can  affect  the  structural  response.  The  third  issue  is  the 
efficient  wall  clock  time  for  simulation,  which  has  to  mainly  rely  on  parallel  computation. 

1.2  The  Objectives 

The  objective  of  this  research  is  to  develop  a  high  accuracy  and  high  efficiency  tool  to 
simulate  the  non-linearity  of  FSI.  It  is  achieved  by  accomplishing  the  following  tasks: 

1)  Develop  a  DES  strategy  using  high  order  scheme  for  turbulence. 

2)  Develop  an  high  order  WENO  scheme  for  shock  capturing  and  high  order  central 
differencing  schemes  for  viscous  terms. 

3)  Develop  a  high  efficiency  low  diffusion  E-CUSP  scheme  to  be  used  as  the  Riemann 
solver. 

4)  Develop  a  high  scalability  general  parallel  computation  methodology  for  structured 
grid. 

5)  Employ  a  fully  coupled  fluid-structural  interaction  model  with  implicit  time  march¬ 
ing. 


1.3  Review  of  the  State  of  the  Art 

1.3.1  Fluid-Structural  Interaction  Nonlinearity 

Many  efforts  have  been  made  by  using  reduced  order  models(e.g.  frequency  domain  meth¬ 
ods),  time  accurate  nonlinear  CFD(e.g.  time  domain  Navier-Stokes  equations),  and  CSD 
(computational  structured  dynamics)  solutions  to  study  the  nonlinear  fluid- structural  inter¬ 
action  problems.  The  reduced  order  models  (ROM)  are  mostly  based  on  linearized  flow 
equations  and  linear  or  non-linear  structural  models.  The  advantage  of  ROM  is  their  CPU 
efficiency  by  significantly  reducing  the  CFD  model  size.  The  ROM  methods  include  the 
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modal  analysis  of  flow  field  with  proper  orthogonal  decomposition  (POD)  by  Hall  and 
Dowell,  et  al.  [4-7],  Lucia,  Beran  and  Silva  [8,9]  or  the  Volterra  theory  based  method  used 
by  Silva  and  Raveh  [10]. 

The  reduced  order  models  have  made  significant  progress  and  are  used  to  study  many 
problems  in  both  external  airframe  flows  and  internal  turbomachinery  flows(Silva  et  al. 
2004,Sarkar  et  al.  2004,  Beran  et  al.  2004,  Tang  et  al.  2004,Tran  et  al.  2004,  Moffatt  et 
al.2004)  [11-16].  A  typical  conclusion  is  that  the  onset  of  flutter  and  LCO  can  be  gener¬ 
ally  predicted  quite  well  with  the  ROMs,  but  the  vibration  amplitude  due  to  primarily  the 
aerodynamic  nonlinearity  such  as  LCO  is  often  significantly  over-predicted.  Resolving  the 
non-linearity  of  the  structural  response  does  not  help  much  as  indicated  by  Tang  et  al.  [2]. 

Other  than  the  ROMs,  the  other  path  that  is  vigorously  under  development  to  pursue  cal¬ 
culation  of  fluid-structural  interactions  is  to  solve  the  time  accurate  Euler  or  Navier-Stokes 
equations  in  time  domain  with  loosely  or  fully  coupled  linear  or  non-linear  structural  mod¬ 
els.  The  loosely  coupled  model  means  that  the  structural  response  lags  behind  the  flow  field 
solution.  Within  a  time  step  for  the  loosely  coupled  method,  the  structure  solver  calculates 
the  response  after  the  flow  solver  is  converged.  This  kind  of  methods  may  be  limited  to 
first-order  temporal  accuracy  only  regardless  of  the  temporal  accuracy  of  the  individual 
solvers  [17].  The  fully  coupled  model  is  that  the  flow  field  and  structure  always  respond 
simultaneously  by  exchanging  the  aerodynamic  forcing  and  structural  displacement  within 
each  inner  iteration  of  a  time  step.  Obviously,  only  the  fully  coupled  model  is  rigorous  in 
physical  sense. 

Bendiksen  et  al.  [18]  pioneered  the  research  by  using  an  explicit  CFD  code  coupled  with 
a  structural  integrator  based  on  the  convolution  integral  to  obtain  the  flutter  boundary  for  a 
NACA  64 A0 10  airfoil.  The  loosely  coupled  model  between  the  fluid  and  structural  solvers 
include  the  work  of  Guruswamy  [19],  Lee-Rausch  et  al.  [20],  Smith  [21],  Vermeersch 
et  al.  [22],  Darracq,  et.  al  [23],  Prananta,  et.  al  [24],  Bohbot  et  al.  [25],  and  Blom  et 
al.  [26].  Alonso  and  Jameson  developed  a  model  which  is  close  to  the  fully  coupled  method 
with  the  structural  displacement  updated  every  several  fluid  solver  iterations  [27].  The 
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implicit  Runge-Kutta  method  with  multigrid  acceleration  is  employed  for  the  flow  solver  in 
Alonso’s  work  [27]  [28].  In  1997-98,  Morton  and  Melville  et  al.  developed  an  implicit  fully 
coupled  fluid  structural  interaction  model,  which  used  the  Beam- Warming  implicit  scheme 
for  the  flow  solver  coupled  with  a  modal  structural  solver  [29]  [30]  [17].  In  2001,  Liu  et  al. 
developed  a  fully  coupled  method  using  Jameson’s  explicit  scheme  with  multigrid  method 
and  a  finite  element  structural  model  [31].  In  2004-2007,  Chen  and  Zha  et  al.  [32-34] 
developed  a  fully  coupled  fluid-structural  interaction  methodology  using  implicit  Gauss- 
Seidel  iteration,  which  does  not  have  the  factorization  error  of  the  Beam-Warming  implicit 
method  and  hence  allows  larger  time  step.  A  new  high  efficiency  low  diffusion  upwind 
scheme,  Zha  E-CUSP  scheme,  is  employed  to  reduce  CPU  time  [33,35-37]. 

However,  the  results  form  the  time  domain  time  accurate  nonlinear  flow  solver  have  not 
yielded  more  optimistic  results.  In  general,  again,  the  weak  divergence,  flutter  divergence 
and  onset  of  LCO  can  be  captured  and  predicted  quite  well,  but  the  amplitude  of  LCO 
is  still  significantly  over  predicted.  These  can  be  seen  in  the  work  by  Bendiksen(1992), 
Morton  and  Beran  (1999),  Weber  et  al.(2001),  and  Tang  at  al.  [2,3,38,39].  Even  though 
the  inviscid  Euler  solver  can  capture  the  onset  of  LCO,  the  viscous  Navier-Stokes  solvers 
predict  the  LCO  amplitude  closer  to  experiments  than  the  Euler  solver.  The  LCO  amplitude 
is  also  dependent  on  the  turbulence  model  used  as  pointed  out  by  Tang  et  al.  (2003)  [2] 
when  the  NLR7301  airfoil  is  calculated  [40,41].  The  non-linear  structural  model  does  not 
really  produce  better  results  of  LCO  as  indicated  by  Gordnier  in  2003  [42].  The  calculation 
with  wind  tunnel  wall  porosity  yields  better  LCO  amplitude,  but  a  significant  difference 
remains  if  the  actual  wind  tunnel  wall  porosity  is  used  (Castro  et  al.,  2001)  [43]. 

According  to  Bendiksen  [3],  both  the  LCO  and  transonic  dip  may  be  caused  by  aerody¬ 
namic  non-linearity  due  to  shock  wave/turbulent  boundary  layer  interaction  and  the  resulted 
separation.  Schewe  and  Dietz  et  al.  in  2003  and  2004  [40,41]  tend  to  support  this  hypoth¬ 
esis  even  though  there  was  no  clear  evidence  in  their  experiment  that  the  flow  separation 
occurs  with  the  very  small  LCO  amplitude.  Schewe  et  al.  [40]  attribute  two  nonlinearity 
mechanisms  to  the  amplitude-limiting  phenomenon:  l)the  oscillating  shock  strength  and 
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the  coupled  pulsations  of  the  marginally  separated  flow  beneath  the  shock;  2)  the  trail¬ 
ing  edge  separation  they  deduced  from  the  significant  increase  in  the  r.m.s.  value  of  the 
pressure  fluctuation  near  trailing  edge. 

Schewe  et  al.  [40]  gave  the  following  important  observations  and  questions:  First,  the 
LCOs  they  captured  have  very  small  relative  amplitudes  of  the  plunge  on  the  order  of 
2/1000  to  3/1000  of  the  chord  and  the  pitching  amplitude  less  than  1°.  They  questioned  if 
the  LCOs  of  such  small  amplitudes  are  the  artifacts  of  the  wind-tunnel  experiment  or  could 
also  occur  in  the  unbounded  flow.  Second,  since  all  the  reported  numerical  simulations  at 
that  time  only  captured  the  much  greater  LCO  amplitudes,  they  questioned  if  it  means  the 
co-existence  of  multiple  LCOs  at  constant  flow  conditions.  They  discovered  the  multiple 
coexisting  LCOs  in  their  experiments.  Third,  they  found  that  the  wall  boundary  layer  tran¬ 
sition  from  laminar  to  turbulent  does  not  have  much  effect  on  LCO.  Fourth,  they  verified 
that  the  wind  tunnel  wall  interference  with  or  without  perforated  test  section  does  not  have 
much  effect  on  LCOs.  Fifth,  they  observed  that  the  transition  from  steady  to  oscillatory 
state  can  be  either  continuous  or  discontinuous. 

The  question  we  want  to  ask  is:  where  is  the  difficulty  for  fluid-structural  interaction 
calculation?  The  answer  appears  to  be  the  aerodynamic  non-linearity.  Bendiksen  pointed 
out  in  1992  [3]  that  both  the  LCO  and  transonic  dip  may  be  caused  by  aerodynamic  non¬ 
linearity  due  to  shock  wave/turbulent  boundary  layer  interaction  and  the  resulted  sepa¬ 
ration.  The  experiments  of  Schewe  and  Dietz  et  al.  in  2003  and  2004  [40,41]  for  the 
NLR7301  supercritical  airfoil  support  these  findings.  Their  experimental  results  indicate 
that  the  flow  separation  at  the  trailing  edge,  and  the  interactions  between  the  shock  and  the 
marginal  region  of  separated  flow  beneath  it,  may  be  responsible  for  limiting  the  amplitude 
of  the  observed  LCO. 

For  large  structural  deflections  and  rotations  in  structure,  Bendiksen  and  Seber  recently 
point  out  that  both  structural  and  aerodynamic  non-linearities  can  have  a  dramatic  qualita¬ 
tive  and  quantitative  effect  on  the  flutter  behavior  [1],  Neglecting  either  the  structural  or 
the  fluid  non-linearities  can  lead  to  completely  erroneous  stability  predictions. 
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In  this  research,  only  small  deflections  and  rotations  in  structure  are  considered.  There¬ 
fore,  the  linear  structural  equations  are  solved  in  the  fluid- structural  interaction  calcu¬ 
lation.  A  fully  coupled  FSI  model  is  used  to  ensure  the  simulation  accuracy  [32-34]. 
The  NLR7301  airfoil  with  detailed  experimental  measurement  of  limited  cycle  oscillation 
(LCO)  is  calculated  in  2D  using  RANS  method  and  in  3D  using  DES  method.  This  research 
is  the  first  effort  to  apply  the  high  order  accuracy  numerical  algorithms  and  DES  to  resolve 
flow  non-linearity  of  flow-structural  interactions.  The  rigorous  algorithm  of  this  research 
appear  to  be  paid  off  with  the  numerical  simulation  matching  the  experiment  excellently 
for  the  first  time.  This  simulation  also  confirms  some  of  the  experimental  observations  and 
answers  some  important  questions.  First,  the  LCOs  with  the  small  relative  amplitude  is 
captured  with  unbounded  flows  in  the  numerical  simulation.  This  means  they  should  not 
be  the  artifacts  of  the  wind-tunnel  experiment  and  most  likely  are  the  factual  phenomenon. 
Second,  the  co-existence  of  multiple  LCOs  at  constant  flow  conditions  is  confirmed  in  our 
simulation.  The  reason  that  other  numerical  simulations  only  capture  the  LCOs  with  large 
amplitudes  may  be  due  to  their  high  numerical  dissipation  that  either  smears  out  the  small 
amplitude  LCO  or  is  only  able  to  resolve  the  large  amplitudes  LCOs.  Third,  the  numerical 
simulation  of  this  research  confirms  that  the  wall  boundary  layer  transition  from  laminar 
to  turbulent  does  not  have  a  large  effect  on  LCOs  at  high  Reynolds  number,  because  our 
simulation  assumes  that  the  boundary  layer  is  fully  turbulent  from  the  airfoil  leading  edge. 
Fourth,  the  simulation  confirms  that  the  wind  tunnel  wall  interference  with  or  without  per¬ 
forated  test  section  does  not  have  much  effect  on  LCOs  because  our  simulation  uses  the 
unbounded  flow  condition  with  no  wind  tunnel  wall  effect  at  all.  Fifth,  the  numerically 
captured  LCO  is  not  accompanied  with  any  flow  separation  due  to  the  very  small  ampli¬ 
tude.  This  may  rectify  the  hypothesis  that  the  LCOs  are  caused  by  the  nonlinearity  of  flow 
separation  induced  by  shock/boundary  layer  interaction.  In  other  words,  the  nonlinearity 
of  shock/boundary  layer  interaction  with  no  flow  separation  is  sufficient  to  trigger  a  LCO. 
This  may  make  reduced  numerical  models  feasible  to  capture  LCOs. 
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1.3.2  Detached-Eddy  Simulation 

To  resolve  the  aerodynamic  non-linearity  such  as  shock  wave/turbulent  boundary  layer  in¬ 
teraction  and  flow  separation,  the  turbulence  simulation  is  critical.  The  widely  used  meth¬ 
ods  for  predictions  of  turbulent  flow  are  based  on  Reynolds  averaged  Navier-Stokes  equa¬ 
tions  (RANS).  However,  RANS  models  are  not  able  to  calculate  the  flow  separation  cor¬ 
rectly.  RANS  methods  intend  to  calculate  the  large  scale  eddies  using  a  universal  model. 
Large  scale  turbulence  is  affected  by  the  flow  geometry  and  boundary  conditions  and  a 
universal  model  does  not  exist. 

Large  Eddy  Simulation  (LES)  is  promising  to  overcome  the  disadvantages  of  the  RANS 
model.  In  LES,  the  governing  equations  are  spatially  filtered  on  the  scale  of  the  numerical 
grid.  The  large  energy  containing  scales  are  directly  simulated,  and  the  small  scale  eddies, 
which  are  generally  more  homogeneous  and  universal,  are  modeled.  The  large  eddies  are 
strongly  affected  by  the  flow  field  geometry  boundaries.  Therefore  the  direct  computation 
of  the  large  eddies  by  LES  is  more  accurate  than  modeling  the  large  eddies  by  RANS. 
The  effect  of  the  unresolved  small  scales  of  motion  is  modeled  by  a  subgrid-scale  (SGS) 
model  [44]  [45]  [46]  [47]  [48]  or  by  the  inherent  dissipation  in  the  numerical  schemes  [49] 
[50].  Because  the  statistics  of  the  small  scale  turbulence  are  more  isotropic  and  universal, 
a  general  physical  model  for  small  scale  eddies  is  more  plausible. 

However,  for  high  Reynolds  number  flows  such  as  those  of  transonic  wings  and  turbo¬ 
machinery  blades,  to  resolve  the  wall  boundary  layer,  LES  needs  the  CPU  resource  not 
much  less  than  the  Direct  Numerical  Simulation(DNS).  This  makes  the  LES  too  expensive 
for  high  Reynolds  number  flow  calculations.  For  engineering  applications,  it  is  not  hopeful 
for  LES  to  be  rigorously  used  until  in  another  4  decades  [51]. 

To  overcome  the  intensive  CPU  requirement  for  LES,  Spalart  et  al.  developed  the  so 
called  detached-eddy  simulation  (DES)  strategy  [51],  which  is  a  hybrid  RANS  and  LES 
method.  Near  the  solid  surface  within  the  wall  boundary  layer,  the  unsteady  RANS  model 
is  realized.  Away  from  the  wall  surface,  the  model  automatically  converts  to  LES.  By  using 
the  RANS  model  near  the  wall,  the  mesh  size  as  well  as  the  CPU  time  can  be  tremendously 
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reduced.  The  motivation  of  DES  is  that  the  LES  is  powerful  in  regions  of  massive  separa¬ 
tion  and  other  free  shear  flows  such  as  jets,  but  much  too  costly  in  the  large  area  of  the  wall 
boundary  layers. 

Spalart  gave  a  grid  guidance  for  DES  in  2001  [52,53],  which  divides  a  flow  domain 
with  solid  walls  to  Euler  region,  RANS  region,  and  LES  region.  In  the  RANS  region,  the 
domain  is  further  divided  to  Viscous  region  and  Outer  region.  The  LES  region  is  composed 
of  Viscous,  Focus  and  Departure  region.  Spalart’s  DES  grid  guidance  gives  sufficient  grid 
resolution  for  LES  region  and  the  transition  to  Euler  region  from  RANS  region.  The  grid 
size  is  dramatically  reduced  compared  to  the  pure  LES. 

Even  though  DES  concept  is  much  newer  than  RANS  and  LES  concept,  its  application 
for  turbulence  simulation  has  already  achieved  encouraging  success  as  shown  in  the  work 
of  Tarvin  et  al.  (1999)  [54],  Spalart  (2001)  [52,53],  Forsythe  et  al.(2002)  [55],  Viswanathan 
et  al.  [56],  Squires  et  al.  [57,58],  Hsnsen,  et  al.  (2003)  [59],  Subbareddy  et  al.  (2005)  [60], 
Wang  et  al.  (2008,  2009)  [61,62].  These  flows  calculated  using  DES  include  those  for 
airfoils,  cylinders,  forbodies,  base  flows,  etc.  The  results  are  qualitatively  and  quantitatively 
better  than  the  solutions  using  RANS.  DES  appears  to  be  a  suitable  compromise  between 
the  physical  models  of  turbulence  and  CPU  efficiency.  In  those  DES  applications,  almost 
all  the  algorithms  use  2nd  order  accuracy  except  that  of  Tarvin  et  al.  (1999)  [54],  which 
employs  fifth  order  upwind  scheme  for  the  inviscid  convective  terms  in  space,  but  2nd  order 
accuracy  for  the  viscous  terms. 

However,  with  the  spread  of  the  successful  DES  application  after  it  was  first  proposed  in 
1997,  a  defect  of  the  first  generation  DES  model(DES97)  [51]  has  been  also  exposed.  It  is 
that  the  transition  from  the  RANS  model  to  LES  in  DES  97  is  not  grid  spacing  independent 
[63].  DES  is  originally  designed  to  treat  the  entire  boundary  layer  using  a  RANS  model 
and  to  use  LES  separated  flow  regions.  A  fine  mesh  with  grid  spacing  much  smaller  than 
the  boundary  layer  thickness  may  exhibit  an  incorrect  behavior  in  boundary  layers  and 
shallow  separation  regions  due  to  locating  the  RANS/LES  transition  within  the  boundary 
layer.  The  grid  spacing  could  be  fine  enough  for  the  DES  length  scale  to  follow  the  LES 
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branch,  which  will  lower  the  eddy  viscosity  below  the  RANS  level.  The  resolved  Reynolds 
stresses  determined  from  the  velocity  fluctuation  (LES  content)  may  be  lacking  because 
the  resolution  is  not  fine  enough  to  fully  support  it.  The  DES  limiter  then  reduces  the 
eddy  viscosity,  and  therefore  the  modeled  Reynolds  stresses.  This  phenomenon  is  referred 
as  modeled-stress  depletion  (MSD)  [63].  This  drawback  is  also  considered  as  one  of  the 
possible  causes  for  the  inaccurate  prediction  of  flow  separation  region  size  with  suction 
flow  control  when  the  DES  is  used  as  indicated  by  Rumsey  [64]. 

To  overcome  the  MSD  problem  and  make  the  DES  limiter  independent  of  grid  spac¬ 
ing,  Spalart  suggested  a  modification  to  the  original  DES97  model  in  2006  [63],  referred 
to  as  Delayed  DES(DDES).  A  blending  function  similar  to  the  one  used  by  Menter  and 
Kuntz  [65]  for  the  SST  model  is  introduced  to  limit  the  DES  length  scale  to  ensure  the 
transition  of  RANS  to  LES  be  independent  of  grid  spacing.  The  DDES  model  has  demon¬ 
strated  excellent  agreement  with  experiment  and  a  significant  improvement  over  the  DES  97 
for  the  tested  cases,  which  include  a  flat  plat  boundary  layer  resolved  with  mesh  spacing 
significantly  smaller  than  the  boundary  layer  thickness,  a  circular  cylinder,  a  single  airfoil 
with  weak  separation  near  trailing  edge,  the  backward  facing  step  with  large  separation 
region,  and  a  multi-element  airfoil.  The  predicted  separation  onset  and  separation  region 
length  agree  well  with  the  experiments. 

This  research  has  successfully  implemented  the  Spalart’s  DES  model,  which  is  used  for 
a  turbulent  cylinder  flow  [61],  prediction  of  NLR7301  LCO  and  CFJ  airfoil  flows  at  high 
angle  of  attack. 

Due  to  the  time  limitation,  this  research  only  adopts  the  DES97  model  [51].  The  DDES 
will  be  implemented  in  future. 

1.3.3  High  Order  WENO  Scheme 

Developing  accurate  and  efficient  numerical  schemes  is  one  of  the  most  important  tasks 
of  the  CFD  researchers  and  engineers.  It  is  particularly  important  when  a  high  fidelity 
numerical  simulation,  such  as  detached-eddy  simulation  (DES)  and  large  eddy  simulation 
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(LES)  is  performed  for  a  unsteady  flow  problem,  which  is  usually  very  CPU  intensive. 
So  far,  most  of  the  engineering  applications  employ  the  2nd  order  numerical  accuracy. 
The  high  order  schemes  (higher  than  3rd  order)  are  mostly  limited  to  the  fundamental 
research  such  as  high  fidelity  turbulence  simulation  (e.g.  Large  Eddy  Simulations  and 
Direct  Numerical  Simulation)  and  aeroacoustic  calculation.  The  reason  is  that  the  high 
order  schemes  are  generally  not  mature  enough  for  robust  engineering  applications.  For 
example,  when  a  high  order  scheme  is  used,  the  convergence  for  a  steady  state  solution  is 
usually  difficult.  How  to  make  a  high  order  scheme  converge  well  is  not  well  studied. 

For  aerospace  engineering  applications  with  shock  waves  or  contact  surfaces,  the  essen¬ 
tially  non-oscillatory  (ENO)  or  weighted  essentially  non-oscillatory  (WENO)  schemes  are 
attractive  for  their  capability  to  capture  discontinuities  and  achieve  the  consistent  high  or¬ 
der  accuracy  in  smooth  regions.  By  using  a  convex  combination  of  all  candidate  stencils  to 
replace  the  smoothest  one  in  the  ENO  scheme,  a  WENO  scheme  has  more  advantages  over 
its  ENO  counterpart.  For  example,  it  approaches  certain  high  order  accuracy  in  smooth 
regions  and  has  better  convergence  rate  due  to  the  smoother  numerical  flux  used.  From 
its  appearance  [66,  67]  to  present,  the  WENO  schemes  have  been  extensively  applied  to 
different  flow  problems  in  many  areas. 

Titarev  and  Toro  [68]  firstly  carried  out  an  extension  of  the  finite-volume  WENO  schemes 
to  three  space  dimensions  with  high  order  accuracy.  A  finite  volume  WENO  scheme  needs 
higher  computational  cost  than  a  WENO  finite  differencing  scheme.  A  WENO  finite  differ¬ 
ence  method  is  more  efficient  in  multi-dimensional  calculation  due  to  avoiding  the  Gaus¬ 
sian  integrals.  As  pointed  out  in  references  [68,69],  when  the  piece-wise  parabolic  re¬ 
construction  is  used  in  two  space  dimensions,  a  finite  volume  WENO  scheme  require  ap¬ 
proximately  three  times  more  CPU  time  than  the  corresponding  finite  difference  WENO 
schemes.  In  three  space  dimensions,  the  difference  is  about  nine  times.  Hence,  for  struc¬ 
tured  meshes,  the  finite  difference  WENO  scheme  is  preferred.  In  [70-72],  the  formally 
high-order  accurate  WENO  shock-capturing  method  using  a  third-order  total- variation  di¬ 
minishing  (TVD)  Runge-Kutta  time  evolution  scheme  is  applied  to  the  re-shocked  two- 
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dimensional  single-mode  Richtmyer- Meshkov  instability  [70],  the  shallow  water  and  the 
open-channel  flow  equations  [71],  and  to  study  adaptive  mesh  refinement  techniques  for 
multi-dimensional  hydrodynamic  simulation  [72].  Sjogreen  and  Yee  [73]  used  a  low  dis¬ 
sipation  sixth-order  spatial  and  fifth-order  WENO  scheme  with  the  standard  fourth-order 
Runge-Kutta  method  to  study  the  supersonic  reactive  flows. 

In  a  WENO  scheme,  a  Riemann  solver  is  needed  to  capture  the  discontinuities.  There 
are  two  ways  to  evaluate  the  Riemann  solver  fluxes.  For  WENO  finite  difference  schemes, 
Shu  suggested  that  the  WENO  reconstruction  be  directly  applied  to  the  split  fluxes  from 
left  or  right  [74].  In  this  research,  we  employ  a  different  method,  which  is  to  evaluate 
the  conservative  variables  with  WENO  scheme  and  then  use  the  conservative  variables  to 
calculate  the  fluxes  based  on  the  Riemann  solvers.  This  is  similar  to  the  MUSCL  method 
of  Van  Leer  [75]. 

Chen  et  al.  [76]  presented  a  class  of  implicit  WENO  schemes  for  the  incompressible 
Navier-Stokes  equations,  in  which  the  lower-upper  symmetric  Gauss-Seidel  (LU-SGS) 
relaxation  is  used  for  computing  steady  state  solutions.  Yang  et  al.  [77]  have  extended 
this  method  to  the  solutions  of  steady  compressible  Navier-Stokes  equations.  Cadiou  and 
Tenaud  [78]  proposed  an  implicit  WENO  shock  capturing  scheme  for  unsteady  flows  and 
applied  it  to  one-dimensional  Euler  equations.  The  use  of  WENO  spatial  operator  not  only 
enhances  the  accuracy  of  solutions,  but  also  improves  the  convergence  rate  of  the  steady 
state  computation  compared  with  using  the  ENO  counterpart.  In  references  [79,  80],  it  is 
found  that  the  factored  LU-SGS  is  significantly  less  efficient  than  the  unfactored  Gauss- 
Seidel  line  relaxation  method  for  the  steady  state  flow  computation  since  the  former  intro¬ 
duces  the  factorization  error  limiting  the  CFL  number  and  convergence  rate. 

Zhang  and  Shu  [81]  found  that,  when  a  5th  order  WENO  scheme  combined  with  a 
Runge-Kutta  time  discretization  is  used  to  achieve  steady  state  solutions,  the  residual  stops 
dropping  at  the  truncation  error  level  of  the  scheme,  which  is  far  above  the  machine  zero. 
They  noticed  that  the  original  smoothness  indicator  of  Jiang  and  Shu  [67]  results  in  a  small 
oscillation  near  a  steady  shock  wave.  The  oscillation  propagates  to  the  smooth  region  and 
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causes  the  residual  to  hang  at  the  truncation  error  level  rather  than  to  approach  machine 
zero.  They  proposed  a  modified  smoothness  indicator  near  the  shock  region  for  the  fifth 
order  WENO  scheme,  which  can  drive  the  residual  to  machine  zero  for  some  ID  and  2D 
problems  without  the  influence  from  the  boundary  conditions.  But  for  the  other  examples, 
the  residuals  still  fluctuate  at  the  level  of  (10  2  ~  10  4).  Zhang  and  Shu  [81]  attribute  the 
convergence  difficulty  to  the  influence  of  boundary  conditions.  At  a  critical  point  (the  first 
derivative  is  zero),  the  first  term  in  the  Taylor  series  expansion  of  the  IS^  of  Zhang  and  Shu 
does  not  satisfy  the  requirement  of  ISk  =  D(  1  +  0(Axr  '))  to  achieve  5th  order  accuracy. 
Thus  the  accuracy  of  the  scheme  of  Zhang  and  Shu  [81]  is  only  3rd  order  at  a  critical  point. 

Henrick  et  al.  [82]  proposed  a  mapped  WENO  scheme  to  achieve  the  optimal  accuracy 
order  at  the  critical  point  of  a  smooth  function  and  discussed  the  choice  of  £  value  for  the 
5th  order  WENO  scheme.  When  £  is  dominant  in  magnitude,  the  preconditions  of  WEN05 
scheme  approaches  those  of  a  central  difference  scheme.  Furthermore,  the  oscillation  on 
the  order  of  £2  may  exist  near  discontinuities.  Hence  if  the  £  is  too  large,  it  will  mitigate 
the  ENO  behavior  of  the  method.  Henrick  et  al.  suggested  that  £  can  be  slightly  larger  than 
the  square  root  of  the  smallest  positive  number  allowed  for  a  particular  machine.  But  they 
didn’t  study  the  convergence  behavior  for  computing  steady  state  solutions. 

So  far,  the  5th  order  WENO  schemes  are  mostly  used  for  unsteady  flow  calculation  such 
as  LES,  DNS,  or  aeroacoustic  calculation  [83-85].  For  unsteady  calculation,  if  an  explicit 
scheme  such  as  a  Runge-Kutta  scheme  is  used,  the  convergence  is  generally  not  an  issue.  If 
a  dual  time  stepping  procedure  [86]  is  used,  the  convergence  could  be  an  issue  within  each 
physical  time  step.  However,  for  dual  time  stepping,  a  fixed  number  of  iterations  within 
each  physical  time  step  is  often  used  and  the  convergence  of  the  solution  is  sometimes 
overlooked.  The  best  convergence  test  is  to  calculate  steady  state  solutions  to  see  if  they 
can  be  converged  to  machine  zero.  For  the  transonic  flows  with  shock  wave  discontinuities, 
little  research  has  been  done  to  study  the  convergence  behavior  for  steady  state  solutions 
using  the  high  order  WENO  schemes. 

In  this  research,  the  5th  order  finite  differencing  WENO  scheme  [87]  is  used  to  evaluate 
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the  inviscid  fluxes,  and  the  4th  order  central  differencing  scheme  [88]  is  used  to  calculate 
the  viscous  fluxes.  The  WENO  scheme  adopted  uses  an  optimized  e  value  and  is  able 
to  remove  the  weights  oscillation,  maintain  the  sensitivity  to  shock  and  contact  surface 
discontinuities,  achieve  optimal  weights  and  thus  the  minimal  dissipation,  and  obtain  solid 
convergence  to  machine  zero. 

1.3.4  Upwind  Schemes 

An  upwind  scheme  is  required  as  a  Riemann  solver  when  a  WENO  scheme  is  used.  Dif¬ 
ferent  from  the  central  differencing  schemes  [89],  upwind  schemes  are  designed  to  make 
the  flux  computation  based  on  the  characteristic  directions  of  the  governing  equations.  The 
upwind  schemes  have  inherent  numerical  dissipation,  which  makes  the  artificial  dissipation 
unnecessary.  However,  if  the  numerical  dissipations  is  too  large,  the  physical  dissipation 
can  be  distorted  [90] . 

The  approximate  Riemann  solver  scheme  developed  by  Roe  [91]  is  one  of  the  most 
famous  upwind  schemes.  By  introducing  the  Jacobian  and  Roe’s  average  for  the  variables, 
the  Roe  scheme  exactly  satisfies  the  Rankine-Hugonoit  relations  and  directly  capture  the 
discontinuities.  The  Roe  scheme  was  considered  as  the  most  accurate  scheme  among  the 
available  differencing  schemes  in  1987  [92].  But  the  Roe  scheme  uses  matrix  dissipation 
and  hence  it  is  time  consuming. 

To  achieve  the  purpose  of  efficiency,  accuracy  and  simplicity  to  use,  many  efforts  have 
been  made  to  develop  upwind  schemes  only  using  scalar  dissipation  instead  of  matrix  dissi¬ 
pation  such  as  that  of  the  Roe’s  flux  difference  splitting  (FDS)  scheme  [91].  The  examples 
include  AUSM  family  schemes  of  Liou  [93-97],  the  Van  Leer-Hanel  scheme  [98],  Ed¬ 
wards’s  LDFSS  schemes  [99, 100],  Jameson’s  CUSP  schemes  and  limiters  [101-103],  and 
the  E-CUSP  schemes  developed  by  Zha,  et  al.  [36,37, 104-108],  etc. 

Pioneered  by  Liou  and  Steffen  [93,95,96],  the  researchers  seeking  the  scalar  dissipation 
primarily  follow  the  guideline  that  the  velocity  and  pressure  should  be  separated  to  consider 
their  characteristics  representing  the  physics  of  the  convection  and  waves.  Liou  and  his 
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colleagues  termed  their  schemes  as  advection  upstream  splitting  method(AUSM)  schemes, 
and  Jameson  gave  the  name  of  convective  upwind  and  split  pressure  (CUSP)  schemes  [101— 
103], 

The  CUSP  schemes  can  be  basically  categorized  to  two  types,  the  H-CUSP  and  E- 
CUSP  [101-103].  The  H-CUSP  schemes  have  the  total  enthalpy  from  the  energy  equation 
in  their  convective  vector,  while  the  E-CUSP  schemes  use  the  total  energy  in  the  convective 
vector.  The  Liou’s  AUSM  family  schemes,  Van  Leer-Hanel  scheme  [98],  and  Edwards’s 
LDFSS  schemes  [99, 100]  belong  to  the  H-CUSP  group.  The  schemes  developed  by  Zha, 
et  al.  [36, 104-106]  belong  to  the  E-CUSP  group. 

Even  though  the  H-CUSP  schemes  such  as  AUSM  family  schemes  have  achieved  great 
success,  from  the  characteristic  theory  point  of  view,  the  schemes  are  not  fully  consistent 
with  the  disturbance  propagation  directions  [37, 109],  which  may  affect  the  stability  and 
robustness  of  the  schemes.  By  splitting  the  eigenvalues  of  the  Jacobians  to  convection 
(velocity)  and  waves  (speed  of  sound),  one  will  find  that  the  convection  terms  only  contain 
the  total  energy  [104],  which  will  lead  to  the  E-CUSP  schemes. 

Zha  and  Hu  recently  suggested  an  E-CUSP  schemes,  which  has  low  diffusion  and  can 
capture  crisp  shock  wave  profiles  and  exact  contact  discontinuities  [36].  The  scheme  is 
consistent  with  the  characteristic  directions  due  to  the  nature  of  E-CUSP  scheme.  The 
scheme  shows  the  highest  stability  for  two  shock  tube  tests  problems  compared  with  several 
other  popularly  used  upwind  schemes  for  the  explicit  Euler  time  marching  scheme.  The 
scheme  also  works  well  when  extended  to  multi-dimensions  [36].  However,  the  E-CUSP 
scheme  of  Zha-Hu  may  generate  temperature  oscillation  near  the  computation  boundary, 
in  particular  when  the  mesh  is  skewed.  Zha  was  able  to  remove  the  temperature  oscillation 
by  introducing  the  total  enthalpy  in  the  smooth  factor  for  the  energy  equation  [37, 109]. 
However,  the  scheme  loses  the  capability  to  capture  the  exact  contact  surface  due  to  the 
modification. 

In  this  research,  a  LDE  scheme  is  developed  by  modifying  the  Zha-Hu  E-CUSP  scheme 
using  the  Mach  number  splitting  of  Edwards’s  LDFSS  schemes  [99, 100]  for  the  convective 
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flux.  The  solutions  calculated  by  the  new  scheme  is  smooth  and  the  scheme  can  capture 
crisp  shock  profile  and  exact  contact  discontinuity. 

1.3.5  Implicit  Time  Marching  Methods 

The  implicit  methods  for  compressible  viscous  flow  calculation  have  been  widely  employed 
due  to  their  less  stiffness  and  faster  convergence  rate  than  the  explicit  schemes.  In  general, 
implicit  methods  require  the  inversion  of  a  linearized  system  of  equations.  The  direct  inver¬ 
sion  of  the  linear  equations  is  usually  preventively  expensive.  The  implicit  linear  equations 
are  therefore  commonly  inverted  by  iterative  methods. 

It  is  known  that  the  approximately  factored  (AF)  implicit  schemes  such  as  the  Beam- 
Warming  scheme  [110]  will  introduce  the  factorization  errors,  which  limit  the  size  of  the 
allowable  time  steps.  For  3D  linear  wave  equation,  the  AF  scheme  is  also  not  uncon¬ 
ditionally  stable.  The  unfactored  schemes  with  no  factorization  errors  such  as  the  line 
Gauss-Seidel  iterations  can  have  larger  time  steps  with  faster  convergence  rate  than  the  AF 
methods  [33,36, 111-114].  However,  the  unfactored  schemes  typically  require  more  CPU 
time  per  iteration  since  the  matrices  are  usually  the  full  Jacobian  matrices  and  can  not  be 
diagonalized. 

The  lower-upper  symmetric  Gauss-Seidel  (LU-SGS)  method  suggested  by  Jameson  and 
Yoon  [115,116]  has  been  widely  used  due  to  their  relatively  easier  implicit  implementation 
[117-119].  The  attractive  feature  of  the  LU-SGS  is  that  the  evaluation  and  storage  of 
the  Jacobian  matrices  can  be  eliminated  by  making  some  approximations  to  the  implicit 
operator.  Although  the  LU-SGS  method  could  be  more  efficient  than  explicit  schemes  and 
is  unconditionally  stable  for  linear  wave  equation,  the  factorization  is  approximated  and 
will  necessarily  introduce  the  factorization  errors. 

For  the  unfactored  implicit  Gauss-Seidel  relaxation  scheme  used  to  solve  the  2D  incom¬ 
pressible  Navier-Stokes  equations,  Rogers  [120]  compared  the  efficiency  of  point-Jacobi 
relaxation  (PR),  Gauss-Seidel  line  relaxation  (GSLR),  incomplete  lower-upper  decomposi¬ 
tion,  and  the  generalized  minimum  residual  method  preconditioned  with  each  of  the  three 
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other  schemes.  If  a  forward  sweep  plus  a  backward  sweep  is  counted  as  one  sweep,  Rogers 
found  that  these  methods  can  obtain  different  efficiency  when  the  different  number  of  the 
sweeps  are  used.  For  three-dimensional  incompressible  flows,  Yuan  [79]  compared  the  ef¬ 
ficiency  of  the  point-Jacobi  relaxation,  line  Gauss-Seidel  relaxation,  and  diagonalized  ADI 
schemes.  Yuan  [79]  observed  that  the  PR(2)  (PR  with  two  sweeps)  is  optimum  in  all  PR(n), 
and  GSLR(l)  is  optimum  in  all  GSLR(n).  For  the  line  Gauss-Seidel  relaxation  methods, 
one  can  choose  one  or  more  of  the  coordinate  directions  as  the  sweep  direction  [1 12, 121]. 
For  compressible  flows,  there  is  few  study  on  how  the  sweep  directions  will  affect  the 
convergence  rate  and  CPU  time.  There  is  also  no  study  to  compare  the  efficiency  of  the 
unfactored  GSLR  and  the  factored  LU-SGS. 

For  unsteady  flow,  Jameson  formulated  a  so  called  dual  time  stepping  scheme  [86].  By 
introducing  a  pseudo  time  term,  the  unsteady  problem  at  each  physical  time  step  is  treated 
as  a  steady  state  problem  for  pseudo  time.  Without  losing  time  accuracy,  the  dual  time 
stepping  scheme  can  greatly  improve  the  computation  efficiency  by  enhancing  diagonal 
dominance  [28].  It  has  been  widely  used  by  researchers  [122-125]. 

In  this  research,  the  unfactored  implicit  Gauss-Seidel  line  relaxation  scheme  is  used  to 
solve  the  3D  compressible  Navier-Stokes  equations.  For  unsteady  computation,  Jameson’s 
dual  time  stepping  scheme  [86]  is  employed  to  facilitate  the  implicit  iteration.  The  temporal 
term  is  discretized  using  the  2nd  order  backward  differencing.  A  comparison  indicates  that 
the  Gauss-Seidel  relaxation  has  higher  convergence  rate  and  CPU  efficiency  than  the  LU- 
SGS  and  GSLR  methods  [80]. 

1.3.6  Parallel  Computation 

Parallel  computing  is  becoming  a  powerful  tool  for  computational  fluid  dynamics  (CFD) 
simulations  [126,  127].  By  interconnecting  PCs  or  workstations  one  can  develop  a  dis¬ 
tributed  parallel  computing  system  to  increase  the  computing  power.  Hence  there  are  many 
efforts  to  develop  parallel  computation  codes  or  to  convert  legacy  sequential  codes  for 
multi-processor  parallel  computation. 
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Multi-block  structured  grids  are  widely  used  for  CFD  parallel  computation  because  of 
their  numerical  efficiency  and  accuracy.  The  basic  idea  is  to  partition  a  large  domain  into 
multiple  smaller  sub-domains  and  conduct  the  computation  in  the  sub-domains  simulta¬ 
neously  to  save  wall  clock  time.  Partitioning  to  multiple  sub-domains  also  decreases  the 
difficulties  of  grid  generation  around  complex  configurations  since  structured  grids  can 
be  generated  within  each  sub-domain  independently.  The  CFD  solvers  written  for  par¬ 
allel  computation  are  usually  based  on  a  single  program  multiple  data  (SPMD)  strategy, 
which  uses  the  same  CFD  code  for  each  sub-domain.  The  sub-domain  data  exchange  at 
domain  partitioning  boundaries,  or  inner  boundaries,  is  usually  treated  as  boundary  condi¬ 
tions  [128].  The  data  are  exchanged  across  the  boundaries  by  a  mapping  procedure  after 
each  iteration  [129].  The  mapping  procedure  determines  the  parallel  computation  effi¬ 
ciency,  accuracy  and  robustness. 

Various  mapping  procedures  for  multi-block  structured  grids  have  been  developed  since 
parallel  computation  was  introduced  to  CFD.  However,  the  procedure  and  implementation 
are  usually  ad-hoc  and  different  code  developers  may  use  different  methods.  The  complex¬ 
ity  of  the  procedures  depends  on  the  complexity  of  the  geometry  to  be  dealt  with.  For  a 
simple  domain  partitioning  problem,  Evans  et  al.  [130]  developed  a  toolkit  known  as  Com¬ 
puter  Aided  Parallelization  Tools(CAPTools)  to  convert  a  scalar  code  to  a  form  suitable  for 
parallel  implementation  with  message  passing  calls.  For  complex  geometries  that  have  dif¬ 
ferent  local  mesh  index  systems,  it  is  difficult  and  time  consuming  to  treat  the  inner  bound¬ 
ary  data  exchange.  One  method  to  avoid  this  is  to  use  database  systems  to  manage  CFD 
parallel  computation.  The  NSMB  package  is  a  well  known  CFD  solver  developed  based  on 
a  database  system  called  MEM-COM  [131, 132].  The  portable  parallel  library  APPL  and 
a  database  manager  GPAR  are  used  to  implement  the  parallel  computation  [133].  How¬ 
ever,  as  database  systems  are  generally  dependent  on  computer  platforms,  portability  of  the 
CFD  codes  are  quite  limited.  In  addition,  such  database  systems  are  not  often  available  for 
general  CFD  code  developers.  Therefore,  it  is  useful  and  necessary  to  develop  a  general 
mapping  procedure  for  inner  boundary  data  exchange  when  a  new  or  legacy  structured  grid 


18 


CFD  code  is  to  be  implemented  for  parallel  computation. 

In  this  research,  a  general  subdomain  boundary  mapping  procedure  is  developed  to 
implement  parallel  computation  for  structured  grid.  A  high  scalability  is  achieved  [134] 

1.4  The  Strategy  of  This  Research 

To  achieve  the  research  objectives,  a  numerical  strategy  to  apply  the  detached-eddy  simu¬ 
lation  (DES)  with  a  fully  coupled  fluid-structural  interaction  model  is  established  for  the 
first  time.  The  following  novel  numerical  algorithms  are  also  created:  a  general  sub-domain 
boundary  mapping  procedure  for  parallel  computation  to  reduce  wall  clock  simulation  time, 
an  efficient  and  low  diffusion  E-CUSP  (LDE)  scheme  used  as  a  Riemann  solver  to  resolve 
discontinuities  with  minimal  numerical  dissipation,  and  an  implicit  high  order  accuracy 
weighted  essentially  non-oscillatory  (WENO)  scheme  to  capture  shock  waves. 

The  Detached-Eddy  Simulation  is  based  on  the  model  proposed  by  Spalart  in  1997. 
Near  solid  walls  within  wall  boundary  layers,  the  Reynolds  averaged  Navier-Stokes  (RANS) 
equations  are  solved.  Outside  of  the  wall  boundary  layers,  the  3D  filtered  compressible 
Navier-Stokes  equations  are  solved  based  on  large  eddy  simulation(LES).  The  Spalart- 
Allmaras  one  equation  turbulence  model  is  solved  to  provide  the  Reynolds  stresses  in  the 
RANS  region  and  provide  the  subgrid  scale  stresses  in  the  LES  region. 

An  improved  5th  order  finite  differencing  weighted  essentially  non-oscillatory  (WENO) 
scheme  with  an  optimized  e  value  is  employed  for  the  inviscid  fluxes.  The  new  LDE 
scheme  used  with  the  WENO  scheme  is  able  to  capture  crisp  shock  profiles  and  exact 
contact  surfaces.  A  set  of  fully  conservative  4th  order  finite  central  differencing  schemes 
are  used  for  the  viscous  terms. 

The  3D  Navier-Stokes  equations  are  discretized  based  on  a  conservative  finite  differ¬ 
encing,  which  is  implemented  by  shifting  the  solution  points  half  grid  interval  in  each 
direction  on  the  computational  domain.  The  solution  points  are  hence  located  in  the  center 
of  the  grid  cells  in  the  computational  domain  (not  physical  domain).  This  makes  it  possible 
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to  use  the  same  code  structure  as  a  2nd  order  finite  volume  method.  A  finite  differencing 
high  order  WENO  scheme  is  used  since  a  finite  differencing  WENO  scheme  is  much  more 
efficient  than  a  finite  volume  WENO  scheme. 

The  unfactored  line  Gauss-Seidel  relaxation  iteration  is  employed  for  time  marching. 
For  the  time  accurate  unsteady  simulation,  the  temporal  terms  are  discretized  using  the 
2nd  order  accuracy  backward  differencing.  A  pseudo  temporal  term  is  introduced  for  the 
unsteady  calculation  following  Jameson’s  method.  Within  each  physical  time  step,  the 
solution  is  iterated  until  converged  based  on  pseudo  time  step. 

A  general  sub-domain  boundary  mapping  procedure  is  developed  for  arbitrary  topol¬ 
ogy  multi -block  structured  grids  with  grid  points  matched  on  sub-domain  boundaries.  The 
interface  of  two  adjacent  blocks  is  uniquely  defined  according  to  each  local  mesh  index 
system  (MIS)  which  is  specified  independently.  A  pack/unpack  procedure  based  on  the 
definition  of  the  interface  is  developed  to  exchange  the  data  in  a  ID  array  to  minimize  data 
communication.  A  secure  send/receive  procedure  is  employed  to  remove  the  possibility  of 
blocked  communication  and  achieve  optimum  parallel  computation  efficiency.  Two  terms, 
“Order’  and  “Orientation” ,  are  introduced  as  the  logics  defining  the  relationship  of  adja¬ 
cent  blocks.  The  domain  partitioning  treatment  of  the  implicit  matrices  is  to  simply  discard 
the  comer  matrices  so  that  the  implicit  Gauss-Seidel  iteration  can  be  implemented  within 
each  subdomain.  This  general  sub-domain  boundary  mapping  procedure  is  demonstrated 
to  have  high  scalability. 

The  validation  computations  are  conducted  for  the  developed  LDE  scheme  and  paral¬ 
lel  algorithm.  These  validation  cases  include  subsonic  flat  plate  turbulent  boundary  layer 
flow,  RAE2822  transonic  airfoil  turbulent  transonic  flow,  transonic  inlet-diffuser  shock 
wave/turbulent  boundary  layer  interaction,  ONERA  M6  Wing  transonic  flow,  3D  Transonic 
Channel  Flow,  Co-flow  jet  (CFJ)  airfoil  internal  and  external  flow  and  circular  cylinder 
massive  turbulent  flow. 

Finally,  the  flow  non-linearity  of  the  fluid- structural  interaction  is  studied  using  Reynolds 
averaged  Navier-Stokes  (RANS)  and  Detached  Eddy  Simulation  (DES)  method  [62].  The 
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fully  coupled  fluid- structural  interaction  procedure  [32-34]  is  employed  to  simulate  the 
FSI  non-linearity  of  NLR7301  airfoil  LCO.  When  using  DES  for  a  2D  geometry,  the  2D 
configuration  is  extended  in  spanwise  direction  to  form  a  3D  geometric  model  since  DES  is 
always  used  in  3D.  The  structural  equations  solved  for  DES  are  the  same  as  those  of  the  2D 
cases.  Two  type  movements  are  included  in  the  present  fluid-structural  interaction:  forced 
pitching  movement  and  flow  induced  vibration. 


Chapter  2 


Governing  Equations 


2.1  3D  Navier-Stokes  Equations 

The  governing  equations  are  the  spatially  filtered  compressible  Navier-Stokes  equations. 
The  spatial  filtering  removes  the  small  scale  high  frequency  components  of  the  fluid  motion, 
while  keeping  the  unsteadiness  associated  with  the  large  scale  turbulent  motion.  Following 
the  derivation  of  Knight  et  al.  [135],  the  filtered  compressible  Navier-Stokes  equations  in 
Cartesian  coordinates  can  be  expressed  as: 

dQ  dE  dF  dG  1  <9EV  dFv  dGv, 

3^  +  3- +  3- +  -5-  = -=-(^  +  -5-^ +  -5-^)  (2.1) 

dt  dx  dy  dz  Re  dx  dy  dz 

where  t  is  time,  Re  is  the  Reynolds  number.  The  variable  vector  Q,  inviscid  flux  vectors  E, 

F,  G,  and  the  viscous  fluxes  Ev,  Fv,  Gv  are  given  as  the  following. 
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The  overbar  denotes  a  regular  filtered  variable,  and  the  tilde  is  used  to  denote  the  Favre 
filtered  variable.  In  above  equations,  p  is  the  density,  u,v,w  are  the  Cartesian  velocity 
components  in  x.y.z  directions,  p  is  the  static  pressure,  and  e  is  the  total  energy  per  unit 
mass. 

The  T  is  the  molecular  viscous  stress  tensor  and  is  estimated  as: 


2  _  dUk  „  _  ,dUj  dUj 

Sii  +  K^r1+  1 
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dxj  dxj 7  ’ 


i,j=  1,2,3 


(2.2) 


The  above  equation  is  in  tensor  form,  where  the  subscript  1,  2,  3  represent  the  coordi¬ 
nates,  x,y,z,  and  the  Einstein  summation  convention  is  used. 

The  molecular  viscosity  jl  =  p  (  f  )  is  determined  by  Sutherland  law. 

The  a  is  the  subgrid  scale  stress  tensor  due  to  the  filtering  process  and  is  expressed  as: 


Gij  =  -P  ( UjUj  -  UiUj) 
The  energy  flux  Q  is  expressed  as: 


(2.3) 


Qi  —  3“  Gij)  Cji  T  (2.4) 

where  <4>  is  the  subscale  heat  flux: 

®i  =  —Cpp(ujT  -  u it)  (2.5) 


The  qi  is  the  molecular  heat  flux: 
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qi  = 


Cpp.  dT 
Pr  dxi 


(2.6) 


pg  =  +l,P(u2+v2  +  w2)  +  pk  (2.7) 

where  7  is  the  ratio  of  specific  heats,  pk  is  the  subscale  kinetic  energy  per  unit  volume. 

pk  =  '-p(wi  -  um)  =  ~\°ii  (2-8) 

In  the  present  calculation,  the  pk  in  Eq.(2.7)  is  omitted  based  on  the  assumption  that 
the  effect  is  small. 

In  generalized  coordinates,  Eq.(2.1)  can  be  expressed  as  the  following: 

SQ  se  av  3&  =  j_/dE1  3K  dG[\ 

»  J5  J?  j;  Re  V  it  Sn  at; ) 

where 


Q=y  (2.10) 

E'  =  j(^Q  +  ^E  +  ^F  +  ^G)  (2.11) 

F;  =  — (t]rQ  +  t\x  E  +  rjy¥  +  r],G)  (2.12) 

G'  =  j(^Q  +  CvE  +  CvF  +  CzG)  (2.13) 

E;  =  i(^Ev  +  ^vFv  +  ^Gv) 


(2.14) 
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F'v  =  j(rix  Ev  +  rjyFy  +  rj- Gv)  (2.15) 

G'v  =  +  C,FV  +  Cz  Gv)  (2.16) 

where  J  is  the  transformation  Jacobian.  The  inviscid  fluxes  in  generalized  coordinate  sys¬ 
tem  are  expressed  as: 
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where  U,  V  and  W  are  the  contravariant  velocities  in  C-  i]  and  £  directions. 
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1,  m,  n  are  the  normal  vectors  on  < § ,  £  surfaces  with  their  magnitudes  equal  to  the  ele¬ 

mental  surface  area  and  pointing  to  the  directions  of  increasing  rj,  £. 
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For  simplicity,  all  the  overbar  and  tilde  in  above  equations  will  be  dropped  in  the  rest  of 
this  report.  Please  note  that  the  Navier-Stokes  equations,  Eq.(2.9),  are  normalized  based  on 
a  set  of  reference  parameters.  The  detailed  normalization  procedure  can  be  found  in  [136]. 
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2.2  Detached  Eddy  Simulation  of  Spalart  [52, 53] 

The  closure  of  the  sub-grid  scale  stresses  and  heat  flux  are  done  by  employing  the  DES  of 
Spalart  [52,53].  The  sub-grid  scale  stresses  are  computed  by: 


. dtij  dUj 
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(2.20) 


The  turbulent  heat  flux  will  be  evaluated  as: 
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where,  pDES  =  pvt. 

The  transport  equation  of  the  Spalart-Allmaras  one  equation  turbulence  model  is  de¬ 
rived  by  using  empiricism,  dimensional  analysis,  Galilean  invariance  and  selected  depen¬ 
dence  on  the  molecular  viscosity  [137].  The  working  variable  v  is  related  to  the  eddy 
viscosity  vf.  The  transport  equation  is  expressed  as 


(2.22) 


757  =  cblSv  (1  -  ft2)  -  [cw\.fw  ~  c^ffa]  [% ] 2 

+  ^[V •  ((v  +  v)Vv)  +  cM(Vv)2]  +  ftl  (Aq)2 

In  generalized  coordinate  system,  the  dimensionlessed  conservative  form  of  Eq.(2.22) 
is  given  as  the  following: 
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The  eddy  viscosity  Vt  is  obtained  from: 

3 

Vt  =  V/V1  fv  1  -  -^-3-  Z  =  l  (2.25) 

X3+<%i  v 

where  v  is  the  kinematic  viscosity.  The  production  term  is: 

s  =  s+^_u,  <2'26) 

where  S  is  the  magnitude  of  the  vorticity.  The  function  fw  is  given  by 

U=g{  )1/6’  S  =  r  +  cw2{rb  -r),  r=-^—  (2.27) 

£6  +  c°3  Skd 

The  function  ft2  is  given  by 

ft 2  =  ct3exp  (-c?4z2)  (2.28) 

and  the  trip  function  f,  \  is 

fn  =  cngtexp  {d2  +  gj dj )  ,  gt  =  min  ^0.1,  (2-29) 

where,  cot  is  the  wall  vorticity  at  the  wall  boundary  layer  trip  location,  d  is  the  distance  to 
the  closest  wall.  dt  is  the  distance  of  the  field  point  to  the  trip  location,  A q  is  the  difference 
of  the  velocities  between  the  field  point  and  the  trip  location,  Avr  is  the  grid  spacing  along 
the  wall  at  the  trip  location. 

The  values  of  the  coefficients  are:  c^i  —  0. 1 355,  cp2  =  0.622,  a  =  \.cw\  =  +  (1  + 

cbl)/(y,cw2  =  0.3,Ch;3  =  2,k  =  0.41, cvi  =7.l,cti  =  l.0,ct2  =  2.0,c?3  =  1.1,Q4  =  2.0. 

In  S-A  one  equation  turbulence  model,  the  trip  point  need  to  be  specified  before  com¬ 
putation.  This  is  not  straightforward  to  do  because  the  exact  position  of  the  trip  point  is  not 
known  in  most  of  the  cases.  Thus,  a  full  turbulent  boundary  layer  is  used  by  setting  ct\  =  0 
and  ct 3  =  0.  No  trip  point  needs  to  be  specified. 

It  is  observed  that  the  S-A  one  equation  turbulence  model  is  sensitive  to  initial  field. 
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If  the  initial  field  of  v  is  set  to  a  small  value,  e.g.  v  <  1,  the  solution  may  converge  with 
v  =  0,  which  is  the  trivial  solution  of  v  when  ct\  —  ct 3  =  0.  This  will  result  in  a  laminar 
flow  solution.  If  the  initial  value  is  too  large  (v  >  3),  the  computation  may  diverge.  In 
addition,  setting  up  the  initial  value  of  v  also  depends  on  the  schemes  to  be  used.  In  our 
computation,  it  is  found  that  it  is  generally  safe  to  set  the  initial  value  of  v  to  2. 

The  boundary  conditions  of  v  are  given  as  the  following 

at  walls  :  v  =  0 

far  field  inflow  :  V  =  0.02 

far  field  outflow  :  V  is  extrapolated 

In  DES,  ct\  =0,  ct 3  =  0.  The  distance  to  the  nearest  wall,  d,  is  replaced  by  d  as 


d  =  min(d,CuES^)  (2.30) 

where  A  is  the  largest  spacing  of  the  grid  cell  in  all  the  directions. 

Within  the  boundary  layer  close  to  walls,  d  —  d,  hence  the  turbulence  is  simulated  by  the 
RANS  mode  determined  by  the  Spalart-Allmaras  model  [137].  Away  from  the  boundary 
layer,  d  =  Cues A  is  most  of  the  cases.  When  the  production  and  destruction  terms  of  the 
model  are  balanced,  the  length  scale  d  will  have  a  Smagorinsky-like  eddy  viscosity  and 
the  turbulence  is  simulated  by  the  LES  model.  Analogous  to  the  classical  LES  theory,  the 
length  scale  A  is  to  cascade  the  energy  to  the  grid  size.  The  coefficient  Cues  =  0.65  is  used 
as  set  in  the  homogeneous  turbulence  [138].  The  Prt  may  take  the  value  of  0.9  within  the 
boundary  layer  for  RANS  mode  and  0.5  for  LES  mode  away  from  the  wall  surface. 

Coupled  Eqs.(2.9)  with  the  S-A  model  Eq.(2.23),  the  conservative  form  of  the  govern¬ 
ing  equations  are  given  as  the  following: 


dQ  dE  d F  dG  _  1  fd R  d S  dT  \ 
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where, 


28 


1 

e  =  7 


P 

pu 

pv 

pw 

pe 

pv 


(2.32) 


pu 

pv 

pW 

puU  +  lxp 

puV  +  mxp 

puW  +  nxp 

pvU  +lyP 

,F  = 

pvV  +  myp 

,G  = 

pvW  +  nyp 

pwU  +  lzp 

pwV  +  mzp 

pwW  +  nzp 

(pe  +  p)U  -ltp 

(pe  +  p)V  -mtp 

(pe  +  p)W  -ntp 

pvU 

pvV 

pvW 

(2.33) 


0 

0 

0 

Ik'Cxk 

nik^xk 

nk^xk 

U'tyk 

,s  = 

mk'tyk 

,T  = 

rik^yk 

U^z.k 

HlkT'zk 

^k^zk 

hfik 

mkPk 

nkPk 

g(v  +  v)(l«Vv) 

^(v  +  v)(m«Vv) 

> 

• 

+ 

Cs.|t> 

_ 1 

(2.34) 
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1 

D  =  - 

J 


0 

0 

0 

0 

0 

Sv 


where,  U,  V,  W  are  defined  as  in  Eq.(2.17). 


(2.35) 


A k  Tki  (Jk  (2.36) 

The  shear-stress  T 'n  and  total  heat  flux  qk  in  Cartesian  Coordinate  can  be  expressed  as 


T ik  =  (H+Ht) 


2  c  dui 
^  Oik  — 

3  dxj 


(2.37) 


7  A-  —  Cp 


\Pr  PrtJ 


dT 

dxk 


(2.38) 


2.3  Structural  Model 

2.3.1  Vortex-induced  Oscillating  Cylinder 

The  structural  model  of  the  vortex- induced  oscillating  cylinder  is  shown  in  Fig.  2.1.  The 
cylinder  is  elastically  supported  and  only  translational  movement  is  considered  in  this 
model. 

The  structural  dynamic  equations  that  govern  the  motion  of  the  cylinder  are: 


mx  +  Cxx  +  Kxx  —  Df 


(2.39) 


my  +  Cyy  +  Kyy  =  Lf 


(2.40) 
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Figure  2.1:  Sketch  of  the  elastically  mounted  cylinder 

These  equations  are  solved  implicitly  together  with  the  equations  of  flow  motion  in  a 
fully  coupled  manner.  In  Eq.  (2.39),  x,  x,  and  x  represent  the  dimensionless  horizontal 
acceleration,  velocity  and  displacement  of  the  moving  object  respectively.  Similarly,  y,  y, 
and  y  in  Eq.  (2.40)  represent  the  acceleration,  velocity  and  displacement  in  the  vertical 
direction.  The  terms  m,  Ly,  and  Df  are  the  mass,  lift,  and  drag  per  unit  span  respectively, 
Cx  and  Cy  are  the  damping  coefficients  in  horizontal  and  vertical  directions,  and  Kx  and  Ky 
are  the  spring  constants  in  horizontal  and  vertical  directions.  In  the  present  study,  this  ‘self- 
excited  oscillators’  is  assumed  to  have  the  same  response  in  both  directions,  i.e.  Cx  =  Cy 
and  Kx  —  Ky. 

If  the  normalization  procedure  is  applied  to  Eqs.  (2.39)  and  (2.40)  by  using  the  same 
reference  scales  of  those  used  for  the  equations  of  flow  motion,  the  following  nondimen- 
sional  equations  are  obtained: 
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x  +  2£ 


2 


(2.41) 


y  +  2dl)y+(l)~y=^-cl  (2.42) 

\u  J  \u  J  flsK 

where  £  is  the  nondimensional  structural  damping  coefficient  calculated  by  C,  —  Cx,y/  [2  sJmKx,y] , 
u  is  the  reduced  velocity  defined  by  u  =  £/«,/ (boo),  b  is  radius  of  the  cylinder,  00  = 
the  mass  ratio  defined  by  ps  —  m/np^b2,  and  cj  and  c/  are  the  drag  and  lift  coefficients 
respectively.  Then  the  equations  are  transformed  to  a  state  form  and  expressed  by: 


M^  +  K-S  =  Q 
dt 


where 


(2.43) 


S  = 


x 

y 

\i  J 


,  M  —  [I],  K  — 


0  -1 

I)2  2C  (1) 

0  0 


0 


0 


0 

0 

0 

'h2 


0 

0 


(!)  2 C(l)  ) 

\uJ  ~  \ U  /  / 


/ 


,  Q  = 


\ 


jlsKCd 


\  Mv7TC/  / 


2.3.2  Flow-Induced  vibration  of  NLR7301  Airfoil 


A  theoretical  model  of  the  structural  dynamics  of  the  test  setup  configuration  with  the  two 
degrees  of  freedom  is  sketched  in  Fig.  2.2. 

The  non-dimensional  form  of  the  equations  governing  the  motion  of  the  two  degree  of 
freedom  based  on  the  model  in  Fig.  2.2  can  be  written  as 


( 


1 


X(x 


—Xa 

r 2 
'a 


d2q 


+  2 


<5/,  o)h  0 

0  r2a  8  a  OOa 


CO2  0  \  =  _2_ 
o  r2(02  )  71  ^ 


J 

(2.44) 
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Figure  2.2:  Model  of  the  structural  dynamics 

where,  xa  is  the  static  unbalance,  ra  is  radius  of  gyration,  coa  is  uncoupled  circular  pitching 
frequency,  ft)/;  is  uncoupled  circular  heave  frequency,  8a  is  Lehr  pitching  damping,  <5/,  is 
Lehr  heave  damping,  fls  is  mass  ratio,  c/  and  cm  are  lift  and  moment  coefficients  respec¬ 
tively.  q  is  defined  by: 


(2.45) 


where  h  and  a  are  the  plunging  and  pitching  displacements  respectively,  «o  is  the  off- wind 
value  of  a.  The  Eqs.(2.44)  are  transformed  to  a  state  matrix  form  and  expressed  as 


dS 

M—  +  K  -  S  —  Q 
dt 


(2.46) 


where 
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The  detailed  derivation  of  Eq.(2.44)  is  given  in  Appendix  A. 


Chapter  3 


The  Numerical  Method 


In  this  chapter,  the  numerical  methods  used  to  discretize  the  governing  equations  are  intro¬ 
duced.  The  flow  governing  equations  are  discretized  using  finite  difference  method  with 
a  fully  implicit  manner.  The  inviscid  fluxes  are  discretized  using  a  newly  developed  low 
diffusion  E-CUSP  scheme  [107, 108].  The  fifth-order  WENO  scheme  [87]  is  used  to  recon¬ 
struct  the  conservative  variables  at  volume  interfaces.  A  set  of  fully  conservative  fourth- 
order  accurate  finite  central  differencing  schemes  for  the  viscous  terms  is  employed  in  this 
research  [87].  These  central  differencing  schemes  are  constructed  so  that  the  stencil  widths 
are  within  the  one  of  the  WENO  scheme.  The  structure  governing  equations  are  discretized 
and  solved  implicitly  in  the  same  manner  to  be  consistent  with  the  flow  governing  equa¬ 
tions. 


3.1  F inite  Difference  Discretization  U sing  Implicit  Method 

The  3D  Navier-Stokes  equations  (2.31)  are  discretized  based  on  a  conservative  finite  dif¬ 
ferencing,  which  is  implemented  by  shifting  the  solution  points  half  grid  interval  in  each 
direction  on  the  computational  domain.  The  solution  points  are  hence  located  in  the  cen¬ 
troids  of  the  grid  cells  in  the  computational  domain  (not  physical  domain).  This  makes  it 
possible  to  use  the  same  code  structure  of  a  2nd  order  finite  volume  method. 
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Using  the  computational  grid  shown  in  Fig.  (3.1)  as  an  example,  the  derivative  of  a 
flux  is  discretized  by  a  finite  difference  method.  Taking  inviscid  flux  E  as  an  example,  the 
discretized  conservative  form  of  its  derivative  can  be  written  as  the  following 


<?E  _  E/+1/2  —  Ef-i/2 
W~  A? 
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Figure  3.1:  Sketch  of  the  computational  grid 

Since  At,  =  1,  At]  —  1,  A£  =  1  are  used  in  the  generalized  coordinate,  the  governing 
Eqs.(2.31)  can  be  written  as  the  following  implicit  form: 


1 

At 


-e“)  +  (Ei+.-E,_,)"+1  +  (F.+, 

=  (Ri+.-R,-.)”+1  +  (Vi-S,_ 


(3.2) 


where  n  and  n  +1  arc  two  sequential  time  levels,  which  have  a  time  interval  of  At.  The  5th 
order  WENO  scheme  with  an  upwind  scheme  Riemann  solver  is  used  for  reconstructing 
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inviscid  fluxes  E.+  i,  F  .+ 1  and  G^+i.  A  fully  conservative  4th  order  central  differencing 
scheme  is  used  to  evaluate  the  viscous  fluxes  R,  S,T. 


3.2  The  Low  Diffusion  E-CUSP  (LDE)  Scheme 


The  basic  idea  of  the  LDE  scheme  is  to  split  the  inviscid  flux  into  the  convective  flux  Ec 
and  the  pressure  flux  Ep.  In  generalized  coordinate  system,  the  flux  E  can  be  split  as  the 
following: 


E 


Ec+Ep  = 


pU  } 

0  \ 

puU 

IxP 

pvU 

lyP 

+ 

pwU 

kP 

peU 

pU 

pvU  ) 

0  J 

(3.3) 


where,  U  is  the  contravariant  velocity  in  £  direction  and  is  defined  as  the  following 


U  is  defined  as: 


U  -  If  lXU  -\~  lyV  4“  IvW 


(3.4) 


U  =  lxU  +  lyV  +  lzW 


(3.5) 


The  convective  term,  Ec  is  evaluated  by 
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Let 


Ec  =  pU 


1 

1 

u 

U 

V 

V 

=  PUf,  f  = 

w 

w 

e 

e 

V  *  J 

V  *  / 

C  =  c 

(l2X+l2y+& 

(3.6) 


(3.7) 


where  c  —  \JyRT  is  the  speed  of  sound.  Then  the  convective  flux  at  interface  i+\  is 
evaluated  as: 


E<  =  Ci  [pLC+f[  +  pRC-fR\  (3.8) 

2  ^ 

where,  the  subscripts  L  and  R  represent  the  left  and  right  hand  sides  of  the  interface. 


Cl  =  \  ( CL  +  CR ) ,  C+  =  a+  (1  +(3l)Ml  -  pLM+  -M+ 

2  2 

C  =  ccR  (1  +  pR)MR  —  fiRMR  +M\ 


ML=^,MR=  gf 

2  2 

aL)R  =  \  [l  ±  sign  {MLjR)] 

/3 l,r  =  - max  [0, 1  -  int{\ML.R\)\ 

Mf  =  Mi  M~  =  Mi  ClrC«r  '  ,  <*> 

Mi=pL8+MR-pR8-M+ 

M^R  =  ±i(MLjt±  l)2 

5±  =  |  { 1  ±  .vi'gn  [5(Ml  +  M^)]} 


(pc2)« 

(pc2)l 


The  pressure  flux,  Ep  is  evaluated  as  the  following 
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EP ,  1  = 


\P 
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/Y 

&+ply 

&+plz 

U  +  C\ 
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&~plx 

£?-ply 

P-Ph 


\ 


\_P 


\ 


U-C 1 

2 

0 


/ 


(3.9) 


L  '  /  R 

The  contravariant  speed  of  sound  C  in  the  pressure  splitting  of  energy  equation  is  con¬ 
sistent  with  U .  It  is  computed  based  on  C  as  the  following. 


C  =  C-lt 


(3.10) 


The  use  of  U  and  C  instead  of  U  and  C  in  the  pressure  splitting  of  energy  to  take  into 
account  of  the  grid  speed  so  that  the  flux  will  transit  from  subsonic  to  supersonic  smoothly. 
When  the  grid  is  stationary,  lt  —  0,  C  =  C,  U  =  U. 

The  pressure  splitting  coefficient  is: 

^t*=  I(M1,s±l)2(2TMt)  (3.11) 

The  LDE  scheme  can  capture  crisp  shock  profile  and  exact  contact  surface  disconti¬ 
nuities.  Since  the  scheme  uses  scalar  dissipation,  for  DES  with  one  extra  equation,  the 
splitting  is  basically  the  same  as  the  original  scheme.  This  is  an  advantage  over  the  Roe 
scheme,  for  which  the  eigenvectors  need  to  be  derived  when  any  extra  equation  is  added 
to  the  governing  equations.  It  is  also  more  CPU  efficient  than  the  Roe  scheme  due  to  no 
matrix  operation. 
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3.3  The  Fifth-Order  WENO  Scheme  [87] 

The  interface  flux,  E.,  1  =  E(Ql- (Jr),  is  evaluated  by  determining  the  conservative  vari- 

1  '  2 

ables  Qi  and  Qr  using  fifth-order  WENO  scheme  [87].  For  example, 


(2l),-+^  =  &»0<70  +  ^1^1  +  ®2<?2 


where 


<7o  =  \Qi- 2  -  \Qi-  \  +  yf  2/ 
<7t  —  —  g2/-t  +  §2,  +  y2/+i 

<72  =  3C?/  T  g2/+l  g2;+2 


05q  +  . . .  +  OJf-_  1 


(3.12) 


(3.13) 


(3.14) 


ak  —  e+isk ,  k  —  0, . . . ,  r  1 

C0  =  0.1,  Ci  =0.6,  C2  =  0.3 

/So  =  if  (Qi- 2  -  2a_i  +  Qi)2  + 1  (Qi- 2  ~  42,-i  4-  3Q,)2  (3.15) 

151  -  ]§  (2,-i  -  22/  +  2/+i)2  + 1  (2/+i  -  2/-i)2 

152  =  If  (Qi~  22/+i  +  G/+2)2  + 1  (32,  -  42, +1  +  2/+2)2 

where,  £  is  originally  introduced  to  avoid  the  denominator  becoming  zero  and  is  supposed 
to  be  a  very  small  number.  In  [87],  it  is  observed  that  IS^  will  oscillate  if  £  is  small  and 
also  shift  the  weights  away  from  the  optimum  values  in  the  smooth  region.  The  higher 
the  £  values,  the  closer  the  weights  approach  the  optimum  weights,  C*-,  which  will  give 
the  symmetric  evaluation  of  the  interface  flux  with  minimum  numerical  dissipation.  When 
there  are  shocks  in  the  flow  field,  £  can  not  be  too  large  to  maintain  the  sensitivity  to  shocks. 
In  [87],  the  optimized  value  of  £  =  10  2  is  recommended  for  the  transonic  flow  with  shock 


waves. 
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3.4  Finite  Differencing  Discretization  of  Viscous  Terms  [87] 

We  take  the  viscous  flux  derivative  in  ^-direction  as  the  example  to  explain  how  the 
schemes  are  constructed.  To  conservatively  discretize  the  viscous  derivative  term  in  Navier- 
Stokes  equations  Eq.(2.31),  we  have 


dR  _  Ri+i/i  —Ri-i/2 

W~  A? 

To  obtain  4th  order  accuracy,  R  needs  to  be  reconstructed  as 

/+ 1/2 

Rj-1/2  —  52  ai 

I=i- 3/2 


(3.16) 


(3.17) 


where 


a, 


1  _  26  1 
-3/2  -  ~24>  ai- 1/2  -  ^  Oi+ 1/2  -  “24 

1/2  =  [(^xTxx)  +  (fly  Try)  +  (CzTrz)]; 

(Ti 


-1/2 


Tju:,/ 


If  in  Eq.(3.17) 


(&|f)+(nz^)  +  (C,|f)]} 


ran  hp 


A(«i+l/2-«i-l/2)  =«'(!, )  +  0(A?4)  (3.19) 

and  the  4th  order  accuracy  is  achieved  (to  be  proved  later).  It  needs  to  point  out  that  in 
Eq.(3.16),  Rj_  1/2  can  not  be  replaced  by  i/2.  Otherwise,  the  4th  order  accuracy  can  not 
be  achieved  even  though  the  high  order  approximation  of  /?,■_ is  used.  The  4th  order 
accuracy  from  Eq.  (3.16)-(3.19)  is  also  based  on  the  uniform  spacing  =  C. 

In  order  to  achieve  the  highest  order  accuracy  of  Rj  with  I  —  i  —  3/2,i  —  1/2,  i  +  1/2, 
the  approximation  of  each  term  in  Eq.  (3.17)  using  the  same  points  is  given  below: 
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M/  =  £  c/jUf+f, 

I=m 

(3.20) 

du  1  ^  I 

w  I  =  ^LD>Ui^ 

(3.21) 

(3.22) 

where 

|  _  1  y1 

At7^  7 

(3.23) 

By  choosing  different  ranges  for  (m,n),  (r.  ,v),  (p,  q)  and  different  coefficients  C[,  D\.  Cj, 

one  can  obtain  different  order  accuracy  approximation  to  the  viscous  terms.  The  principle 

of  choosing  (m,n),(r,s),(p,q)  is  to  ensure  that  the  approximation  of 

|f  If  in  Eq.(3.16) 

is  a  central  differencing.  For  example,  let  (, m,n )  =  (—2, 1),  (r.s)  —  (  — 

3,2), and  ( p,q )  = 

(—2, 2),  and  they  give 

M/=  £c/ju«+/  +  o(a^4), 

f=m 

(3.24) 

^  j^/Wf+f  +  ^A^5), 

(3.25) 

(3.26) 

where 

drj^J  =  Ap  t  C^'j+/  +  °(Ar?4) 

(3.27) 

the  coefficients  CLdLG  can  be  obtained  by  Taylor’s  series  expansion  and  are  given  in 
Tables  3. 1-3.3.  For  example, 


Mi-3/2  =  i^(5M/-2  +  15/if- !  -  5m  +  Hi+i)  +  0(A%4) 

<  Mf-t/2  =  i^(— Mi-2  +  9/if-i  +9/i/  —  /i/+i) +  0(Ai§4)  (3.28) 

k  Mi+1/2  =  i^(M/-2  —  5/if-i  +  15/if +  5/i/+i)  +  (9(Ai§4) 
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( 


< 


v 


du 

d%  I  '-3/2 

du  | 

5?  li-1/2 

eh/  | 

I  i+l/2 


(  1920 Mi— 3  —  T28Ui— ~  64W,— 1  192mi  —  128M/+1  640M/+2)  “1“  b(Ac,  ) 

Af  (  —  640m'“3  +  5^4m/'-2  —  64m/'-1  +  64w/  —  5MW'+1  +  640Mi+2)  +  0(A <^5  ) 

A^(  — 640M/-3  +  I§8M'— 2  —  192M/-1  ~~  6 HM*  +  T28M'+I  —  l920M'+2)  +  0(AE,~  ) 

(3.29) 


The  other  terms  are  determined  similarly.  For  comparison,  the  terms  used  in  Ref.  [139, 
140]  by  De  Rango  and  Zingg  et  al.  are  given  as  the  following, 


Ah'- 3/2  —  T6  (— Ah'-3  +  9/i/-2  +  9/t,_i  —  Ah)  +  0(A/§4) 

<  Ah— 1/2  =  t^(M/-2  +  9^-1  +  9/1,  -  Ai/'+i)  +  <9(A^4)  (3.30) 

,  Ah+t/2  =  i/>(Ah-t  T  9/1;  +  9/i;+ 1  —  /i;+2 )  +  <9(A/§4) 

||li-3/2  —  24^  (— w/'-3  ~27n;_2  +  27n,-_i  —  w;)  +  0(A(§4) 

<  |f|;-l/2  =  24A?(-M«-2-27M;_1+27M;-M;+1)  +  0(A^4)  (3.31) 

,  If  I  i+l/2  —  —  ' Ui~l  ~  27w;  +  27«;+i  —  W;+2)  +  <9(  A<§4) 

Compare  Eqs.(3.28),(3.29)  and  Eqs.(3.30),(3.31),  it  can  be  seen  that  /!/  in  present  paper 
has  the  same  accuracy  order,  as  that  of  De  Rango  and  Zingg  et  al.,  but  has  small  stencil 

width  (i  —  2,  ■■■  ,i  +  1),  ||| /  has  the  same  stencil  width,  but  obtains  one  accuracy  order 

higher  than  that  in  Ref.  [139, 140]. 


Table  3.1:  The  coefficients  of  Cj 


I 

rj 

C  2 

CLi 

C1 

co 

q 

/  —  3/2 

5/16 

15/16 

-5/16 

1/16 

1—1/2 

-1/16 

9/16 

9/16 

-1/16 

i  + 1  /2 

1/16 

-5/16 

15/16 

5/16 

It  can  be  proved  that  the  scheme  Eq.  (3.16)  is  symmetric  with  respect  to  cell  i.  For 
example,  the  coefficients  of  /l,  _2«;-3,  /i;+2«;+3,  /l;  i  «;_2,  and  /i,+im(+2  can  be  found  as 
(in  the  following  formula,  Cj  and  b\  are  the  coefficients  of  /l,  h/,  ui+i  in  R/  for  Ri+  j/2, 
respectively.  It’s  clear  that  there  are  Cj  —  C\z\  and  b\  —  D\  J ,  a/  =  «/_i ,  /  =  i 1/2,  i  4- 
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Table  3.2:  The  coefficients  of  D\ 


I 

#-3 

D- 2 

^-i 

^0 

^2 

z-  3/2 

71/1920 

-141/128 

69/64 

1/192 

-3/128 

3/640 

z  —  1/2 

-3/640 

25/384 

-75/64 

75/64 

-25/384 

3/640 

z  + 1/2 

-3/640 

3/128 

-1/192 

-69/64 

141/128 

-71/1920 

Table  3.3:  The  coefficients  of  Cf 


cc 

C  2 

cii 

rc 

co 

q 

Cc 

c2 

1/12 

-8/12 

0 

8/12 

-1/12 

l/2,i  +  3/2): 


C, 


i—2,i—3 


=  ~t^na,CI_2D 


3/2 


_  A.2Zl_  +  26.('=lw^3\  +  /=l', 

Lv  24  3  16  1920  '  24  1  16  >  \  640 ^24/  16  U40/J 


46080 


Q+2./+3  =Li+Jl_\l2aiCI2DI3 


1/2' 


=  ( =1 ).  2,  .  j-  +  26  .  (=1  \ .  _3_  ,  t  zT\ .  A  .  ( 
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Q-1,1-2  =  uiCi_1D,_2  - 1£%  aiC'_iD'_2 
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So  we  have  Q_2j!-_ 3  =  Q+2/+3,  Q_i)!-_2  =  Q+1)f+2,  and  so  on.  Hence  the  scheme  Eq. 
(3.16)  is  symmetric  with  respect  to  grid  node  i.  The  symmetry  of  central  differencing  for 
Eq.  (3.16)  satisfies  the  diffusion  property  of  the  viscous  flux. 

Next,  we  prove  that  the  order  of  accuracy  given  by  Eq.(3.19)  is  satisfied.  Take  the  term 
T  —  jidu/dt,  in  Eq.(3.19)  as  the  example, 

In  Ri-in,  at  /  =  i  —  3/2,  based  on  Taylor’s  series  expansion 


7 


:3/2  =1  ni=mc\\iM^Y!l=rD\uM) 

=  Mj-3/2  +AijijJ3i2A^4 +  0(A^5)  f||,-_3/2  +  0(A%~ 

=  Mi- 3 /2  ff  I  /— 3/2  +^/M,^3 /2  If  I  i-3/2^  4  +  <^(  A<§  5 ) 


A/  is  the  coefficient  of  Taylor’s  series  expansion. 

The  corresponding  term  7’+  in  Ri+ j /2  is  at  /  =  i  —  1/2,  and 


TA 


1/2  =ILiMC/Mi+l+/(^ELr^+1+/) 

=  Mi-1/2  +^/M,^i/2^4  +  ^(^5)  ffli-l/2  +  0(A^' 

=  Mi- 1  /2 1|  I  i- 1/2  +  ^/M^4 i /2  If  I »- 1  /2  4  +  ^  5  ) 
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Note  that  A/  =  A/,  hence 

du .  du  * 

Ti- 1/2  _  Ti- 3/2  =  I. i-1/2  “  Mi-3/2 li-3/2  +  ) 

The  other  two  terms  can  be  analyzed  similarly  as  above,  then  Eq.(3.19) 

^rW+i/2  ~Ri- 1/2)  =  *'(&)  +  o(A|4) 

is  proved,  i.e.  the  constructed  schemes  are  formally  4th  order  accuracy. 

3.5  Implicit  Time  Integration 

3.5.1  Implicit  Flow  Solver 

The  time  dependent  governing  equations  are  solved  using  dual  time  stepping  method  sug¬ 
gested  by  Jameson  [86].  To  achieve  high  convergence  rate,  the  implicit  pseudo  time  march¬ 
ing  scheme  is  used  with  the  unfactored  Gauss-Seidel  line  relaxation.  The  physical  temporal 
term  is  discretized  implicitly  using  a  three  point,  backward  differencing  as  the  following 
(The  prime  is  omitted  hereafter  for  simplicity): 

dQ  3(2n+l-40"  +  e'i-1 

dt  2A  t  1  j 

where  n  —  1,  n  and  n  +  1  are  three  sequential  time  levels,  which  have  a  time  interval  of 
At.  The  first-order  Euler  scheme  is  used  to  discretize  the  pseudo  temporal  term  to  enhance 
diagonal  dominance.  The  semi-discretized  equations  of  the  governing  equations  are  finally 
given  as  the  following: 


1  1.5  \  /  c)R  \  w+l  ,m 

At  +  At  J  1  \dQj 


gQn+l,m+l  _  j^n+l,m 


3  Q 


n+l,m 


4  Qn  +  Q 


,n— 1 


2At 


(3.33) 
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where  the  At  is  the  pseudo  time  step,  R  is  the  net  flux  evaluated  on  a  grid  point  using  the 
fifth-order  WENO  scheme. 


3.5.2  Gauss-Seidel  Line  Relaxation  [136] 


To  enhance  diagonal  dominance,  a  first  order  scheme  is  used  for  the  implicit  pseudo  tem¬ 
poral  terms.  Following  the  procedure  in  Hu’s  Ph.D.  report  [136],  the  implicit  discretized 
form  of  Eq.  (3.2)  is  written  as  the  following 


+  A+ A  Ql+l  uk  +A-A0R  u  +  B+AQ?+ 
+B  AQ'l+l 


+  C+A£"j1 


-a  nf1+l  — 


ij,k+ 1 


+  C~AQ n. 


i,j,k- 1 


yi+1 

j+I,k 

RHS" 


(3.34) 


RHS"  is  the  summation  of  all  the  terms  on  the  right  hand  side  (RHS)  of  the  equation. 


RHS'!  =  At 
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+  |Ski 


S"_  \ 
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E"  i  -EM_, 

2  l  2 


+  (f';  i-F';  x 

J'  2  J  2 


+  ig;+1-gj_, 

2  K  2 


j  +  D"  ■  At  (3.35) 


Gauss-Seidel  line  relaxation  is  applied  in  each  direction  (i,  j,  k )  and  is  swept  one  time 
step  forward  and  backward  in  each  direction.  For  example,  the  equation  for  Gauss-Seidel 
relaxation  following  lines  along  direction  i  with  the  index  from  small  to  large  is  written  as: 


B  AQ] 


|M+1 
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+  «A2yi  +  S+Ae” 
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(3.36) 


where 
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RHS"  -A+AQ'; 
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(3.37) 
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3.5.3  Implicit  Structural  Solver 

The  structural  equations  (2.43)  and  (2.46)  are  discretized  and  solved  implicitly  in  each 
physical  time  step  in  a  manner  consistent  with  flow  governing  equations  (3.33): 


1  1.5 

—I+—M  +  K 
At  At 


A  5sn+l,r, 
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m+ 1  _ Qii+l,m+l 
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3Sn+]'m-4Sn  +  Sn-] 
2At 


KS‘ 


n+l,m 


(3.38) 


The  fluid- structural  interaction  is  implemented  in  a  fully  couple  manner  [33,34].  Within 
each  physical  time  step,  the  flow  equations  and  structural  equations  are  solved  iteratively 
until  the  prescribed  convergence  criteria  is  satisfied  for  both  flow  and  structural  solver.  Af¬ 
ter  the  convergence  criteria  is  reached,  the  fluid-structural  interaction  goes  to  next  physical 
time  step. 


3.6  Boundary  Conditions 

To  obtain  a  well  posed  solution  of  a  given  flow  problem  by  solving  the  Navier-Stokes  gov¬ 
erning  equation,  Eq.(2.9),  it  is  necessary  to  define  the  boundary  conditions  for  the  problem. 
Since  the  solution  points  are  located  at  the  centroids  of  the  cells,  the  ghost  cells  are  used  to 
define  the  boundaries  except  for  the  inviscid  flux  on  wall  surface  for  steady  state  problems. 
In  other  words,  most  of  the  boundary  conditions  are  defined  by  setting  up  the  values  of 
the  variables  at  the  ghost  cells.  Depending  on  the  scheme  order  of  accuracy  to  be  used, 
the  number  of  ghost  cells  will  vary  to  match  the  accuracy  of  the  inner  points.  Several 
commonly  used  boundary  conditions  are  described  as  the  following. 

3.6.1  Supersonic  Inflow 

For  supersonic  inflow  boundary,  all  the  primitive  variables  are  fixed  as  the  initial  values  of 
the  flow  field  at  the  ghost  cells. 
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Pgst  —  Pint  ■>  Ugst  —  Mint i  VgSt  —  Vint->  ^gst  —  W int •>  &gst  —  &int  \5.Dy) 

where,  gst  represents  the  ghost  cell  and  int  represents  the  initial  value.  In  this  case,  the 
initial  values  are  set  to  the  values  of  free  stream  and  used  to  specify  the  inflow  BCs. 

3.6.2  Supersonic  Outflow 

For  supersonic  outflow  boundary,  all  the  primitive  variables  are  extrapolated  with  zero 
gradient  from  their  inner  counterparts. 


Pgst  —  Pinn^Ugst  —  ^inn^gst  —  ^inn^gst  —  W  inn-,  &  gst  —  &inn  V-3.4U) 

where,  inn  represents  the  inner  counterpart. 

3.6.3  Subsonic  Inflow  (For  External  Flows) 

For  subsonic  inflow  boundary,  four  variables  are  specified  using  the  free  stream  values  of 
the  flow  field.  One  variable  is  extrapolated  with  zero  gradient  from  its  inner  counterpart. 
Usually,  the  static  pressure  is  extrapolated 


Pgst  —  Pint  i  UgSt  —  Uint ,  VgSf  —  Vjnt ,  WgSt  —  Wjnt ,  Pgst  —  Pinn 


(3.41) 


3.6.4  Subsonic  Outflow  (For  External  Flows) 

For  subsonic  outflow  boundary,  the  static  pressure  is  fixed  as  the  initial  values  of  the  flow 
field.  The  normalized  outflow  static  pressure  is  determined  by  free  stream  Mach  number. 
Pint  —  jtp  ■  F°ur  variables  are  extrapolated  with  zero  gradient  from  their  inner  counterparts. 


Pgst  —  Pinn-t^-gst  —  ^-inni^gst  —  ^inni^gst  —  W inmPgst  —  Pint 


(3.42) 


The  other  variables  are  calculated  based  on  these  5  variables. 
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3.6.5  Subsonic  Inlet  (For  Internal  Flows) 

For  subsonic  inlet  boundary,  the  prescribed  variables  are  usually  total  pressure  pt,  total 
temperature  Tt  and  the  two  flow  angles  a  and  /3.  The  total  pressure  and  temperature  are 
used  because  they  are  easier  to  obtain  in  wind  tunnel  experiments.  A  velocity  component  is 
extrapolated  from  its  inner  counterpart.  The  other  two  velocity  components  are  calculated 
using  the  flow  angles  a  and  /3 .  The  extrapolated  velocity  component  is  determined  by  main 
flow  direction  which  is  given  by  user.  For  example,  if  x  direction  is  taken  as  the  main  flow 
direction,  the  velocity  components  are  calculated  as  the  following: 


Ugst  —  Minn  :  Vgst  —  //(.,/  tan  (x .  WgSt  —  Mgj/tan/3  (3.43) 

The  other  variables  are  calculated  using  the  following  normalized  relations: 


TgSt  —  Tt  —  -  (y-  1)4^,  (ugSt  +  vgst  +WgSt) 


(3.44) 
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Tgst  +  -z  (ugst  +  V2gst  +  wgst) 


y(y-  \)Ml  8sr  2  ™  gst  gst 


(3.47) 


3.6.6  Subsonic  Outlet  (For  Internal  Flows) 

For  subsonic  outlet  boundary,  it  is  the  same  as  the  subsonic  outflow  boundary  except  that 
the  static  pressure  is  fixed  as  a  given  value  poutiet- 


Pgst  —  Pinn ■  Mgst  —  Mjnn.  VgS(  —  V 'lnn .  W gSf  —  Winn:  Pgst  —  Poutiet 


(3.48) 
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3.6.7  Wall  Boundary 

Stationary  Walls 

For  inviscid  flux,  the  wall  boundary  condition  is  enforced  by  setting  the  normal  contravari- 
ant  velocity  on  the  boundary  to  zero.  For  example,  if  a  wall  boundary  is  on  a  ^-surface,  the 
contravariant  velocity  V  is  zero  on  the  wall  surface,  and  hence  the  flux  on  the  wall  surface 
is  calculated  as, 


pV 

0 

puV  +  mxp 

m.xp 

pvV  +  myp 

= 

triy.p 

pwV  +  mzp 

mzp 

(, pe  +  p)V 

W 

0 

The  wall  pressure  is  extrapolated  from  inner  points  by  the  following  formulation: 
1)  1st  order  extrapolation 


2)  3rd  order  extrapolation 


Pw  =  P  t 


(3.50) 


Pw  =  ^(llpi -7p2-2p3)  (3.51) 

o 

For  viscous  flux,  no-slip  and  adiabatic  wall  boundary  condition  is  conducted  by  setting 
the  ghost  cell  velocity  as  the  negative  of  the  velocity  of  its  inner  counterpart  based  on  the 
reflection  condition. 


Pgst  —  Pinn :  Mgst  —  M-inniVgst  —  ^ inn ■,  Wgst  —  Winn,SgSt  —  Sinn 


(3.52) 
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Moving  Walls 

At  2D  moving  boundary  surface,  since  the  wall  velocity  is  not  zero,  both  the  inviscid  and 
viscous  fluxes  on  the  wall  are  evaluated  using  ghost  cells.  The  density  and  velocity  are 
calculated  by  extrapolation. 


Pgst  —  Pinn 

UgSt  —  2.v^  Uj,m  (3.53) 

Vgst  —  2yb  Vjnn 

where,  Xb  and  yt,  are  the  wall  boundary  velocity  components  in  x  and  y  direction  respec¬ 
tively. 

The  inviscid  normal  momentum  equation  is  solved  to  obtain  pressure.  The  momentum 
equation  can  be  described  as  the  following: 


pa—  — Vp  (3.54) 

The  normal  momentum  equation  is: 


Vp-n  =  —  pa  n 


where,  a  is  the  acceleration: 


(3.55) 


a  =  xbi  +  ybj 


(3.56) 


n  is  the  unit  normal  vector: 


n  =  nxi  +  nyj 


1 


(flxi  +  %j) 


(3.57) 


Thus, 
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-pan 


^  ,1/2  (^xXb  +  ^yyh) 
(fll  +  fly) 


(3.58) 
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•  (nj  +  nyj) 
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If  (^  +  ^y)  +  |f  (Vx+Vy) 


<9/? 


(3.59) 


If  the  grid  on  wall  surface  is  orthogonal,  ^Tr]Y  +  <gvr/v  =  0.  Then  Eq.(3.59)  can  be 
expressed  as 


Vp  ■  n  = 


.1/2 


(flX2  +  0y) 

The  normal  momentum  equation  Eq.(3.55)  can  be  written  as  the  following. 


(3.60) 


dp  (172  +  172)  («  + w)  (3-61) 

3.6.8  Symmetrical  Boundary 

For  a  symmetrical  boundary,  the  corresponding  velocity  component  in  the  ghost  cell  is 
mirror  reflected  about  the  symmetrical  boundary  from  its  inner  counterpart.  The  other 
variables  are  extrapolated  with  zero  gradient.  For  example,  if  the  symmetrical  boundary  is 
a  .i'-planc.  the  boundary  condition  is  defined  as  the  following: 


Pgst  —  Pinn ■  Mgst  —  ^inrii  Vgst  —  f inn ■,  W gst  —  ^inn^gst  —  ^inn 


(3.62) 
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3.6.9  Division  Periodic  Boundary 

For  division  periodical  boundary,  the  variables  and  coordinates  are  periodically  defined  in 
the  periodical  direction  (Fig.  3.2).  For  example,  if  the  periodic  boundary  is  in  ^-direction, 
the  variables  are  calculated  in  the  same  way  as  inner  points. 


Figure  3.2:  periodic  boundary 
For  the  ghost  cell  of  the  start  surface  in  ^-direction, 


Pgst  —  PiencljUgst  —  ^iend  i^gst  —  ^iendi^gst  —  ^iendi^gst  —  &iend 

For  the  ghost  cell  of  the  ending  surface  in  ^-direction, 


(3.63) 


Pgst  —  Pistart  i^gst  —  distort  ->V gst  —  V istart  gst 


^  istart-,^  gst  —  &  istart 


(3.64) 


where,  the  istart  and  tend  represent  the  first  and  last  cell  number  in  ^-direction. 
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3.6.10  Translational  Periodic  Boundary  for  Flow  Variables  only 

This  boundary  condition  is  used  for  detached-eddy  simulation  to  simulate  a  2D  isometric 
problem  (Fig.  3.3)  with  spanwise  periodic  BC.  Different  from  the  division  periodic  bound¬ 
ary,  the  geometric  information  of  the  ghost  cells  in  the  translational  periodic  boundary 
are  extrapolated  from  adjacent  inner  points.  The  flow  variables  still  satisfy  Eq.(3.63)  and 
(3.64). 


Y 


Translational  periodic  BC. 


Figure  3.3:  translational  periodic  boundary 

3.6.11  Subdomain  Boundary  (Inner  Boundary) 

This  boundary  condition  is  used  to  define  the  interface  of  two  block  grids  partitioned  for 
parallel  computing  or  multi-block  grid  computation.  The  data  of  the  ghost  cells  or  halo 
cells  are  given  by  exchanging  the  boundary  data  of  the  two  adjacent  blocks.  The  subdomain 
boundary  mapping  procedure  and  data  exchanging  process  are  described  in  Chapter  5. 


Chapter  4 


Comparison  of  the  Low  Diffusion 
E-CUSP  Scheme  with  the  Roe  Scheme 


To  demonstrate  the  accuracy,  efficiency,  and  robustness  of  the  new  LDE  scheme,  several 
2D  and  3D  cases  are  computed  using  the  LDE  scheme  and  the  Roe  scheme  to  compare  their 
performance.  Both  S-A  one  equation  model  and  B-L  algebraic  model  are  used  for  compar¬ 
ison.  The  finite  volume  solver  using  3rd  order  MUSCL  scheme  for  inviscid  fluxes  and  2nd 
order  central  differencing  scheme  for  viscous  terms  is  employed  for  this  comparison  [136]. 

4.1  Subsonic  Flat  Plate  Turbulent  Boundary  Layer  Flow 

The  subsonic  flat  plate  is  used  to  examine  the  performance  of  the  LDE  scheme  for  turbulent 
boundary  layer.  The  mesh  size  is  181  x  81  (see  Fig.  4.1).  The  y+  of  the  first  cell  center  to 
the  wall  is  kept  less  then  1.0.  The  Reynolds  number  is  4  x  106  based  on  the  length  of  the 
plate.  The  inlet  Mach  number  is  0.5. 

As  shown  in  Fig.  4.2,  both  the  computational  results  of  the  LDE  scheme  and  the  Roe 
scheme  agree  well  with  the  law  of  the  wall  using  S-A  model.  They  are  slightly  better  than 
the  results  using  B-L  model  in  the  transition  region  from  the  linear  viscous  sublayer  to  log 
layer.  With  the  B-L  model,  both  schemes  can  use  a  large  CFL  number(>  100).  With  the  S- 
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56 


Figure  4.1:  The  mesh  for  subsonic  flat  plate 

A  model,  the  CFL  number  can  be  set  to  100  for  the  LDE  scheme.  But  for  the  Roe  scheme,  it 
can  be  only  set  to  about  10.  That  means  the  Roe  scheme  needs  more  time  steps  to  converge 
a  result  than  the  LDE  scheme  and  hence,  need  more  computational  time.  Fig.  4.3  shows 
the  solution  residuals  of  the  LDE  scheme  and  the  Roe  scheme  for  S-A  and  B-L  turbulence 
model.  The  LDE  scheme  only  use  about  ^  of  the  iterations  required  by  the  Roe  scheme. 


Figure  4.2:  Comparison  of  velocity  profiles 


Figure  4.3:  The  L2  solution  residual  history  of  flat  plate 
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4.2  RAE2822  Transonic  Airfoil 

The  RAE2822  transonic  airfoil  is  used  to  compare  the  performance  of  the  LDE  and  Roe 
scheme  with  S-A  model  for  computing  2D  turbulent  transonic  flows.  The  mesh  is  a  two- 
block  O-grid  with  dimensions  of  2  x  (129  x  56)  as  shown  in  Fig.  4.4.  The  Reynolds  number 
is  6.5  x  106  based  on  the  chord  length.  The  Mach  number  is  0.729.  The  angle  of  attack  is 
2.31°. 


Figure  4.4:  2-Block  grids  for  RAE2822 

Fig.  4.5  shows  the  convergence  histories  of  the  FDE  scheme  and  the  Roe  scheme.  The 
maximum  CFF  number  that  the  Roe  scheme  can  use  is  6.0,  whereas  the  FDE  scheme  can 
use  10.  The  FDE  scheme  achieves  significantly  faster  convergence  rate  and  lower  residual 
level.  Fig.  4.6  presents  the  comparison  of  pressure  coefficients  between  the  experimental 
data  and  computation  results.  The  results  of  the  FDE  scheme  and  Roe  scheme  are  virtually 
identical  and  the  predicted  shock  locations  agree  well  with  the  experiment. 
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Iterations 


Figure  4.5:  The  L2  solution  residual  history  of  RAE2822 


Figure  4.6:  The  surface  pressure  coefficient  distribution  of  RAE2822 

4.3  Transonic  Inlet-Diffuser 


The  transonic  inlet-diffuser  is  used  to  examine  the  performance  of  the  LDE  scheme  for 
shock  wave/turbulent  boundary  layer  interaction.  The  mesh  is  a  single-block,  2D  H-grid 
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with  dimension  of  193  x  97(Fig.  4.7).  The  Reynolds  number  is  4.38  x  105  based  on  the 
throat  height.  The  inlet  Mach  number  is  0.5.  The  exit  back  pressure  equals  to  0.72  times 
of  the  inlet  total  pressure. 


Figure  4.7:  The  mesh  for  2D  Inlet  diffuser 


Fig.  4.8  presents  the  comparison  of  the  experimental  data  and  the  computational  results. 
It  shows  that,  with  the  S-A  model,  the  LDE  scheme  and  Roe  scheme  have  nearly  identical 
results.  The  S-A  model  predicts  the  results  significantly  better  than  the  B-L  model. 


Figure  4.8:  Upper  wall  pressure  distribution  of  the  inlet  diffuser 

Fig.  4.9  is  the  pressure  contours  of  the  LDE  scheme.  A  curved  A -shock  is  clearly 
captured  due  to  the  shock  wave/turbulent  boundary  layer  interaction. 


Fig.  4. 10  is  the  computed  stream  lines.  It  indicates  that  the  upper  wall  boundary  layer  is 
separated  due  to  the  shock/boundary  layer  interaction.  This  case  shows  that  the  turbulence 


Figure  4.9:  The  contours  of  the  pressure  for  inlet  diffuser 


model  is  a  critical  factor  for  the  prediction  accuracy  of  the  shock  wave/turbulent  boundary 
layer  interaction. 


Figure  4.10:  The  stream  lines  of  inlet  diffuser 


4.4  Transonic  ONERA  M6  Wing 

The  transonic  ONERA  M6  wing  is  calculated  to  examine  the  performance  of  the  LDE 
scheme  for  three  dimensional  cases.  The  mesh  is  composed  of  16  block  grids  which  are 
obtained  by  partitioning  a  single  block  O-H-grid  with  the  dimensions  ofl45x61x41 
(Fig.  4.11).  The  Mach  number  is  0.8395.  The  Reynolds  number  is  1.97  x  107  based  on  the 
averaged  chord.  The  angle  of  attack  is  3.06°. 

Fig.  4.12  and  Fig.  4.13  present  the  comparison  of  the  pressure  distributions  between  the 
experiment  and  computation  at  the  different  sections.  The  location  of  z/b  —  0.2  is  near  the 
root  and  z/b  —  0.99  is  near  the  tip  of  the  wing.  The  computation  results  agree  well  with  the 
experimental  data  except  at  the  section  of  z/b=0.8,  where  the  double-shock  pattern  is  not 
well  resolved  as  most  of  other  CFD  simulations. 

Fig.  4.14  and  Fig.  4.15  plot  the  contours  of  pressure  on  the  surfaces  of  M6  wing  using 
the  FDE  and  Roe  schemes  respectively.  The  pressure  contours  of  the  pressure  and  suction 
surface  are  put  on  the  left  and  right  respectively  to  have  a  clear  view  of  the  3D  wing  sur¬ 
faces.  Both  schemes  clearly  capture  the  flow  pattern  that  two  shock  waves  merge  near  the 
wing  tip  on  the  suction  surface  and  highlight  a  typical  lambda-shape.  One  shock  wave  is 
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Figure  4.1 1:  M6  wing  mesh 

near  the  leading  edge,  and  the  other  impinges  the  root  just  after  the  half  of  the  cord,  while 
it  touches  almost  at  the  leading  edge  at  the  tip.  The  two  figures  indicate  that  the  LDE  and 
Roe  schemes  predict  nearly  the  same  results. 

Similar  to  the  2D  cases,  the  maximum  CFL  number  of  the  LDE  scheme  is  larger  than 
that  of  the  Roe  scheme.  Fig.  4.16  shows  the  convergence  histories  of  the  LDE  scheme  with 
CFL  =  5  and  the  Roe  scheme  with  CFL  =  2.  The  LDE  scheme  achieves  significantly  faster 
convergence  rate  and  lower  residual  level. 
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Figure  4.14:  The  contours  of  surface  pressure  for  M6  wing  using  LDE  scheme 
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Figure  4.15:  The  contours  of  surface  pressure  for  M6  wing  using  Roe  scheme 
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4.5  3D  Transonic  Channel  Flow 

The  3D  transonic  channel  flow  is  used  to  examine  the  performance  of  the  LDE  scheme  and 
the  Roe  scheme  for  3D  shock  wave/turbulent  boundary  layer  interaction  problems  with  the 
S-A  model.  A  single  H-grid  with  the  dimensions  of  90  x  60  x  60  is  used  in  the  computa- 
tion(Fig.  4.17).  The  Reynolds  number  is  106  based  on  the  entrance  height.  For  boundary 
conditions,  the  total  pressure,  total  temperature  and  flow  angle  are  fixed  at  the  inlet  and  the 
static  pressure  at  the  outlet  is  adjusted  to  match  the  position  of  the  shock  wave  obtained  by 
the  experiment. 


Figure  4.17:  Transonic  duct  3D  mesh 

Fig.  4.18  to  Fig.  4.20  present  the  comparison  of  the  Mach  number  contours  between  ex¬ 
periment  and  computation  at  three  spanwise  sections.  The  computed  shock  wave  structures 
of  the  FDE  scheme  and  the  Roe  scheme  agree  well  with  each  other  and  are  similar  to  that  of 
the  experiment.  At  the  two  locations  close  to  the  side  wall  at  z  =  60 mm  and  z  —  90 mm,  both 
the  computations  over  predict  the  size  of  separation  zone.  At  the  mid  section  z  —  10mm, 
the  predicted  separation  size  and  pattern  agree  well  with  the  experiment. 


Chapter  5 


Parallel  Computation 


In  this  chapter,  the  sub-domain  boundary  mapping  procedure  used  for  parallel  computa¬ 
tion  is  described  in  details.  The  3D  Navier-Stokes  solver  using  an  implicit  time  marching 
scheme  with  line  Gauss-Seidel  relaxation  is  parallelized  with  the  SPMD  parallel  strategy. 
This  procedure  is  used  for  both  the  2nd  order  and  high  order  schemes.  Several  2D  and  3D 
cases  are  computed  to  test  the  parallel  computing  efficiency  and  robustness  of  the  parallel 
code. 


5.1  The  Mapping  Procedure 

5.1.1  Inner  Boundary  and  Relationship  Between  Adjacent  Blocks 

The  indices  i.  j,  k  are  used  to  express  the  mesh  index  of  an  arbitrary  sub-domain  block.  For 
any  structured  grid  block  with  the  mesh  dimensions  of  imax  =  n  1 ,  jrnax  —  n2.kmax  =  «3,  we 
can  uniquely  define  an  arbitrary  block,  face  and  edge  using  the  two  diagonal  points  at  the 
opposite  comer  of  each  entity.  For  example,  the  block,  face  and  edge  shown  in  Fig.  5.1  can 
be  defined  as  the  following 
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block  : 


start  =  1,1,  \}end  =  nl,n2,n3 
start  —  77 1, 1,  l, end  =  nl,n2,n3 
start  =  n\ ,  1 , 1 ,  end  —  « 1 , 1 ,  «3 
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/ace  (/  =  n  1 )  : 
edge(i  =  n\.  j  =  1)  : 


(1  ,n2,n3)  (n1,n2,n3) 


block 


(n1,n2,n3) 


(nl.1,1) 


face 


(nl  ,1  ,n3) 

(nl  ,1,1) 
edge 


Figure  5.1:  Definition  of  a  block,  face  and  edge 
Where,  start  and  end  are  one-dimensional  arrays  with  3  elements.  For  simplicity,  as¬ 
sume  that  the  blocks  are  numbered  from  1  to  n  and  there  are  only  two  halo  layers  of  over¬ 
lapping  grid  points  for  the  inner  boundaries  (Fig.  5.2).  Thus  two  layers  of  data  need  to  be 
communicated  at  each  inner  boundary  between  the  adjacent  blocks.  The  number  of  halo 
layers  can  be  arbitrary  depending  on  the  accuracy  order  of  the  scheme  to  be  used.  In  our 
code,  we  have  used  up  to  four  halo  layers  for  7th  order  WENO  scheme  [141]. 

For  an  arbitrary  block  p,  to  uniquely  define  its  inner  boundaries  and  their  relationship 
with  the  adjacent  blocks,  the  following  information  is  needed: 
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(I,n2,n3)  (nl,n2,n3) 


Halo  points  of  block  q 


Cjnterface^ 


(I,n2,n3)  (ml,n2,n3) 


Halo  points  of  block  p 


Figure  5.2:  Inner  boundary  of  two  adjacent  blocks 

1)  The  block  number  of  the  adjacent  block  q. 

2)  The  inner  boundary  definition,  i.e.  the  two  diagonal  points  of  the  boundary  between 
block  p  and  q. 

3)  The  relationship  of  the  MISs  between  block  p  and  q,  which  is  needed  for  pack¬ 
ing/unpacking  the  data  for  exchange. 

A  term,  “ Order ”,  is  introduced  to  express  the  MIS  of  a  block.  The  numbers  1,2,3 
represent  the  mesh  axis  directions  of  i.  j,  k  respectively  for  a  block.  The  mesh  axis  index 
sequence  of  the  current  block  p  will  be  always  (1,2,3).  The  Order  of  the  adjacent  block 
q  is  always  defined  based  on  the  mesh  axis  directions  of  the  current  block  p.  There  are 
6  possibilities  for  the  Order  of  block  q  as  shown  in  Fig.  5.3.  An  Order  is  independent  of 
axis  directions.  For  example,  in  Fig.  5.3,  Order  1  of  the  block  q  shows  two  opposite  I-axis 
direction.  Hence,  for  each  order  there  are  in  total  six  combinations  of  index  axis  directions. 
In  this  manner,  mesh  axis  direction  does  not  need  to  follow  the  right  hand  rule.  The  Order 
is  numbered  from  1  to  6  as  given  in  table  5.1  and  Fig.  5.3.  An  inner  boundary  and  its 
relationship  with  the  adjacent  block  can  be  uniquely  defined  by  the  current  block  number 
p,  the  adjacent  block  number  q ,  the  inner  boundary  diagonal  points  of  block  p  and  block  q, 
and  the  Order  of  block  q. 
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Figure  5.3:  MIS  relationship  of  adjacent  blocks 

The  Order  of  block  p  corresponding  to  the  MIS  of  block  q  can  be  uniquely  defined 
based  on  the  Order  of  q  corresponding  to  block  p.  The  relationship  is  given  in  table  5.2 
based  on  the  relationship  indicated  in  Fig.  5.3.  Table  5.2  indicates  that  most  of  the  orders 
are  the  same  except  Orders  4  and  5. 

With  the  Order  information  of  any  two  adjacent  blocks,  an  inner  boundary  only  needs 
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Table  5.1:  The  Order  of  the  block  q 


Order  Number  of  Block  q 

1 

2 

3 

4 

5 

6 

The  Order  of  Block  q 

(1,2,3) 

(1,3,2) 

(2,1,3) 

(2,3,1) 

(3,1,2) 

(3,2,1) 

Table  5.2:  The  relative  relationship  of  Orders  for  two  blocks 


Order  Number  of  Block  q 

1 

2 

3 

4 

5 

6 

The  Order  of  q  based  on  MIS  of  p 

(1,2,3) 

(1,3,2) 

(2,1,3) 

(2,3,1) 

(3,1,2) 

(3,2,1) 

The  Order  of  p  based  on  MIS  of  q 

(1,2,3) 

(1,3,2) 

(2,1,3) 

(3,1,2) 

(2,3,1) 

(3,2,1) 

Order  Number  of  Block  p 

1 

2 

3 

5 

4 

6 

to  be  defined  for  one  block,  as  the  other  will  take  its  Order  according  to  table  5.2.  Thus, 
the  probability  of  making  mistakes  in  defining  the  inner  boundary  conditions  is  minimized. 
As  an  example,  the  inner  boundary  (see  Fig.  5.2)  between  block  p  and  block  q  can  be 

defined  as  the  following: 

block  —  p, start  =  n  1, 1,  l, end  =  nl,n2,n3, 

iblock  =  q, istcirt  =  1,1,1  pend  —  \,n2,n3,order  —  1,2,3 
Where,  the  block  and  iblock  represent  the  current  block  number  and  adjacent  block 

number  respectively.  The  start,  end  and  istcirt,  iend  are  the  diagonal  points  defining  the 

inner  boundary.  The  start,  end  and  istcirt,  iend  are  given  according  to  the  local  MIS  of  the 

sub-domain,  which  is  independent  of  other  sub-domains.  The  order  in  the  above  example 

represents  the  Order  of  block  q  corresponding  to  block  p. 

5.1.2  Pack  and  Unpack  Data  Procedures 

For  CFD  parallel  computation  with  multi-block  grids  the  inner  boundary  data  of  the  flow 
field  are  exchanged  after  each  iteration  according  to  the  Order  relationship  of  the  inner 
boundary  defined  in  the  last  section.  To  be  efficient,  a  one  dimension  array  is  used  for  data 
communication.  For  each  block,  two  operations  need  to  be  done  for  data  exchange: 

1)  pack  the  inner  boundary  data  into  an  one-dimensional  array  and  send  them  to  the 
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adjacent  block  to  be  used  for  the  halo  cells. 

2)  unpack  the  one-dimensional  array  received  from  the  adjacent  block  and  assign  it  to 
the  halo  cells  of  the  inner  boundary. 

The  pack/unpack  procedure  for  a  2D  problem  is  simpler  than  the  3D  one  as  the  2D 
boundaries  are  edges  and  the  3D  boundaries  are  faces.  Hence,  we  will  describe  the  pack/unpack 
procedures  for  2D  and  3D  problems  separately. 

5.1.2.1  2-D  Problems 

The  pack/unpack  procedures  are  implemented  using  the  following  rules: 

1)  Inner  boundary  data  are  packed  into  an  one-dimensional  array  in  the  inward  direction 
of  the  interface,  i.e.  from  the  outermost  edge  to  the  innermost  edge. 

2)  The  one-dimensional  array  received  from  the  adjacent  block  is  unpacked  in  the  re¬ 
versed  (outward)  direction. 

For  example,  assuming  the  number  of  the  halo  layer  is  l  and  letting  /?3  =  1,  the  pack/unpack 
procedure  of  the  2D  problem  shown  in  Fig.  5.2  for  block  p  is  given  as  the  following: 

Pack: 

il  =0 

do  i  —  start(l),start(l)  —  l  +  1,  —  1 

do  j  —  st art (2),  end (2) 

il  =  il  +  1 

bcb(il)  —  x(i,j ) 

end  do 

end  do 
Unpack: 
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/I  =0 

do  i  =  start(  1)  +  1  .start{\ )  + 1 
do  j  —  st art (2),  end (2) 
il  =  i  1  +  1 
x(i,j)  =  bcb(il ) 
e/7J  Jo 
do 

where,  l  =  2  is  used  for  the  example  shown  in  Fig.  5.2,  ^arf(l)  =  n\, start (2 )  = 
l,e«J(2)  =  /?2,  Z?cZ?  is  the  one-dimensional  array  used  for  data  communication,  and  x  is 
the  data  array  of  the  flow  field. 

5.1.2.2  3-D  Problems 

For  3D  problems  the  data  will  be  packed  from  the  outermost  plane  to  innermost  plane  in 
the  current  block.  The  data  is  unpacked  in  the  reversed  direction.  However,  since  the  inner 
boundary  is  a  surface  for  a  3-D  problem,  the  Orders  of  the  adjacent  blocks  are  needed  to 
determine  the  sequence  that  the  data  is  to  be  packed  and  unpacked  in  the  remaining  two 
directions  on  an  inner  boundary  surface. 

To  define  the  data  packing  sequence,  it  is  necessary  to  introduce  another  term,  “ Orien¬ 
tation ’,  to  specify  the  orientation  and  location  of  the  inner  boundary  interfaces  of  a  block. 
The  Orientation  is  defined  based  on  each  inner  boundary  interface  as  given  in  table  5.3. 
Therefore,  any  inner  boundary  of  a  block  has  an  Orientation  numbered  from  1  to  6. 

Table  5.3:  The  relationship  between  the  Orientation  and  faces 


Face 

i=l 

i=nl 

j=l 

j=n2 

k=l 

k=n3 

Orientation 

1 

2 

3 

4 

5 

6 

Unpacking  in  3D 

First,  we  need  to  decide  a  rule  for  unpacking  for  the  current  block  to  simplify  the 
pack/unpack  procedure.  The  rule  we  adopt  is  that,  when  the  Orientation  is  determined,  the 
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data  unpacking  sequence  for  the  remaining  two  axis  directions  is  always  from  the  larger 
axis  number  to  the  smaller  one.  In  this  way,  the  unpack  procedure  will  be  independent  of 
the  Order  of  the  adjacent  blocks  when  the  Orientation  is  determined.  For  example,  if  the 
Orientation  of  an  inner  boundary  of  the  current  block  p  is  2,  the  unpacking  is  done  by  the 
following  DO  loops: 

Unpack: 

il  =0 

do  i  —  start(l)  +  l,start(l)  + 1 
do  j  —  st art (2),  end (2) 
do  k  =  st  art  (3),  end  (3) 
il  =  il  +  1 
x(i.j.k)  =  bcb(i\) 
end  do 
end  do 
end  do 

If  the  Orientation  of  an  inner  boundary  of  the  current  block  p  is  4,  the  unpacking  is 
done  by  the  following  DO  loops: 

Unpack: 

il  =0 

do  j  =  start( 2)  +  l,  start  (2)  +  1 
do  i  =  start(l),end(l) 
do  k  =  start (3),  end (3) 
il  =  il  +  1 
x(i,j,k )  =  bcb(i  1) 
end  do 
end  do 
end  do 

If  the  Orientation  of  an  inner  boundary  of  the  current  block  p  is  6,  the  unpacking  is 


done  by  the  following  DO  loops: 
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Unpack: 

il  =0 

do  k  =  start (3)  +  1,  start (3)  + 1 
do  i  =  st art (1), end (1) 
do  j  —  st  art  (2),  end  (2) 
il  =  il  +  1 
x(i,j,k)  =  bcb(i  1) 
end  do 
end  do 
end  do 

For  the  odd  Orientation  numbers,  the  DO  loops  of  the  unpack  procedure  are  similar. 

Packing  in  3D 

To  match  the  unpack  procedure,  the  data  packing  sequence  in  the  current  block  must 
match  the  MIS  of  the  adjacent  block  which  is  defined  by  the  Order  of  the  adjacent  block. 

Assuming  that  the  Orientation  of  the  inner  boundary  of  the  current  block  p  is  2,  the 
remaining  two  directions  of  the  data  packing  for  the  current  block  p  are  then  (2,3).  Based 
on  Table  5.1  and  Fig.  5.3,  for  Order  numbers  1,  3  and  5,  the  corresponding  two  directions 
of  the  adjacent  block  <7  are  (2,3),  (1,3)  and  (1,2)  respectively.  To  match  the  data  unpacking 
sequence  in  the  adjacent  block  q ,  the  packing  sequence  in  the  current  block  then  must  be 
from  the  larger  axis  number  to  the  smaller  one.  Specifically,  the  data  packing  DO  loops  are 
the  following: 

Pack: 
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/I  =0 

do  i  =  st  art  {l) ,  start  (l)  —  l  +  1 
do  j  =  st art (2),  end (2) 
do  k  =  st  art  (3),  end  (3) 
il  =  il  +  1 
bcb(i  1)  =  x(i.j.k) 
end  do 
end  do 
end  do 

Where,  i,  j  and  k  correspond  to  1,  2  and  3  respectively.  For  Order  numbers  2,  4  and 
6,  the  corresponding  two  directions  of  the  adjacent  block  q  are  (3,2),  (3,1)  and  (2,1) 
respectively.  To  match  the  data  unpacking  sequence  in  the  adjacent  block  q,  the  packing 
sequence  in  the  current  block  must  be  from  the  smaller  axis  number  to  the  larger  one,  which 
will  have  the  following  DO  loops: 

Pack: 

il  —  0 

do  i  =  st  art  (l) ,  start  (l)  —  l  +  1 
do  k  =  start (3),  end (3) 
do  j  =  start (2), end (2) 
il  =  il  +  1 
bcb{il)  =  x(i.j.k) 
end  do 
end  do 
end  do 

The  above  two  types  of  DO  loops  cover  all  the  scenarios  under  Orientation  2.  For  other 


Orientation  numbers,  the  DO  loops  are  similar. 
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5.1.3  Send/Receive  Procedure 

As  shown  in  Fig.  5.4,  the  send/receive  procedure  for  CFD  parallel  computation  is  for  all 
blocks  to  send  data  to  all  adjacent  blocks  first,  and  then  all  blocks  receive  data  from  all  adja¬ 
cent  blocks  [142].  This  procedure  can  avoid  the  communication  deadlock,  but  a  communi¬ 
cation  blockage  may  occur  because  of  limits  to  computer  buffer  space  which  is  used  to  tem¬ 
porarily  save  the  exchanged  data.  To  avoid  communication  blockage,  a  secure  send/receive 
procedure  is  implemented  in  the  following  way. 

The  basic  idea  of  the  secure  send/receive  procedure  is  to  do  the  send/receive  operations 
in  a  pair  simultaneously.  That  is,  when  a  block  sends  data,  the  receiving  block  will  receive 
the  data  at  the  same  time  (Fig.  5.5).  The  procedure  is  implemented  based  on  the  following 
rules  for  a  block: 

1)  Send  data  to  the  adjacent  blocks  with  greater  block  numbers  and  receive  data  from 
the  adjacent  blocks  with  smaller  block  numbers. 

2)  Send  data  to  the  adjacent  blocks  with  smaller  block  numbers  and  receive  data  from 
the  adjacent  blocks  with  greater  block  numbers. 


Figure  5.4:  The  send/receive  procedure  that 
may  create  buffer  space  blockage 

This  procedure  ensures  the  one-to-one  correspondence  between  send  and  receive  oper¬ 
ations  requiring  minimal  buffer  space  to  avoid  communication  blockage.  Specifically,  the 
procedure  is  implemented  by  two  DO  loops. 
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Figure  5.5:  The  secure  communication  pro¬ 
cedure  that  minimize  buffer  space  usage 

Loop  1  (For  all  interfaces): 

1)  Obtain  the  block  number  of  the  adjacent  block. 

2)  If  the  block  number  of  the  adjacent  block  is  greater  than  the  block  number  of  the 
current  block,  send  data;  otherwise,  receive  data. 

Loop  2  (For  all  interfaces): 

1)  Obtain  the  block  number  of  the  adjacent  block  again. 

2)  If  the  block  number  of  the  adjacent  block  is  smaller  than  the  block  number  of  the 
current  block,  send  data;  otherwise,  receive  data. 

Fig.  5.6  shows  the  flow  chart  of  the  exchanging  procedure. 

This  method  is  proved  to  be  very  effective  in  removing  the  communication  blockage 
problem  and  has  a  high  efficiency  of  data  communication  in  our  numerical  experiments. 

5.2  Implementation 

5.2.1  Flow  Charts 

The  in-house  Navier-Stokes  code  [34-36]  is  converted  to  have  parallel  computing  capabil¬ 
ity  using  the  suggested  general  sub-domain  boundary  mapping  procedure.  Based  on  multi¬ 
block  grids,  the  Single  Program  Multiple  Data  (SPMD)  parallel  strategy  is  employed  in  the 
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Figure  5.6:  Flow-chart  of  exchanging  procedure  for  parallel  computation 

parallelization.  In  this  strategy,  the  CFD  solver  is  designed  as  a  block  solver  which  solves 
flow  governing  equations  in  each  block.  The  interface  is  taken  as  a  boundary,  called  inner 
boundary.  After  a  time  step,  the  data  of  the  halo  cells  of  inner  boundaries  are  exchanged 
across  the  interface  by  the  mapping  procedure  described  in  section  5.1.  The  SPMD  flow 
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chart  for  the  parallel  computation  is  given  as  the  following 


start 


j 


call  mpi_initial 


process  if 


read  datain 
input  parameters 


3 


process 


HE 


read  datain 
input  parameters 


E 


process  rr 


read  datain 


call  initflow 

read  initial  field  and  BC 

call  initflow 

read  initial  field  and  BC 

i  i 

call  exchange 
exchange  halo  cell's  data 

call  exchange 
exchange  halo  cell's  data 

input  parameters 

7 


call  initflow 
read  initial  field  and  BC 


E 


call  integrate_ALL 
iterate  one  time  stepl 


I 


E 


I 


call  exchange 
[exchange  halo  cell's  data) 


call  exchange 
bxchange  halo  cells  datal 


call  integrate_ALL 
iterate  one  time  step 

V 

call  exchange 
exchange  halo  cells  datal 


T 


T 


net 


residual<eps)  (residual<eps 


^no 


call  integrate_ALL 
[iterate  one  time  step 

V 

call  exchange 
bxchange  halo  cells  data) 


3 


no 


residual<eps 


3 


yes 


i 


3 


i 


yes 


call  output 

call  output 

call  output 

output  results 

output  results 

output  results 

3 


call  mpi_final 

call  mpi_final 

call  mpi_final 

Figure  5.7:  The  flow  chart  for  parallel  computation  in  SPMD  strategy 
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5.2.2  Implicit  Gauss-Seidel  Iteration 

Gauss-Seidel  line  iteration  needs  a  special  treatment  when  a  single  block  is  splitted  into 
multi-block  to  do  parallel  computation.  For  example,  following  lines  along  direction  i  with 
the  index  from  small  to  large,  Eq.  (3.36)  can  be  written  in  a  matrix  form  (ignore  the  dash 
lines  here): 


1 

COi 

Co 

1 

AQi 

1 - 

a 
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1 _ 

b2  b2 

Bt  | 
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AQi 
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B,n  Brn  5+ 
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— 
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D  + 

Dm+ 1 

AQm+ 1 

rhs;„+1 

0 

1 

1 

Bn 

Bn 

A  Qn 

RHS' 

When  using  multi-block  parallel  computation,  for  example  two  blocks,  Eq.  (5.1)  will  be 
partitioned  into  two  sub-matrices  as  indicated  by  the  dash  lines  in  Eq.  (5.1).  The  variables 
AQi  — >  A Qm  and  AQm+i  — >  A Qn  are  solved  on  two  separate  processors  by  conducting  the 
matrix  inversion  iteration  on  each  sub-matrix.  A  simple  treatment  is  to  discard  the  comer 
matrices  5+  on  the  first  sub-matrix  and  B~+l  on  the  second  sub-matrix.  They  are  the 
coefficients  computed  based  on  the  variables  from  the  adjacent  sub-domain.  In  the  present 
work,  these  two  coefficients  are  treated  as  zero  in  the  implicit  solver  for  the  sub-domain 
computations.  The  two  sub-domain  matrix  systems  obtained  are: 


Bi  B+ 
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Bm+l  B^n |_1  AQm+l  RHSm+ 1 

Bm+ 2  Bm+ 2  B+ 1_2  AQm+2  _  ^HSm+2 

5-  Bn  AQn  RHSl 

The  advantages  of  discarding  the  matrices  5+  and  B~+l  is  that  it  avoids  exchanging 
the  matrices  across  the  sub-domain  boundaries  and  hence  reduces  communication  time. 
The  disadvantage  is  that  it  does  not  preserve  the  exact  matrices  as  in  the  single  processor 
computation  and  may  affect  the  convergence  efficiency.  However,  the  computation  exper¬ 
iments  indicate  that  the  slowing  down  of  convergence  due  to  this  non-exact  treatment  is 
small,  particularly  when  the  mesh  size  is  large.  It  should  be  pointed  out  that  discarding 
the  matrices  B~+1  and  Bm  will  not  affect  the  accuracy  of  the  solution  when  it  is  converged. 
This  is  because  the  variables  at  the  sub-domain  boundaries  are  exchanged  exactly  match¬ 
ing  the  result  of  the  single  domain  calculation  to  calculate  the  RHS  of  Eqs.  (3.34),  which 
determines  the  accuracy  of  the  solutions. 

5.3  Results  and  Discussion 

To  validate  the  accuracy  of  the  parallel  computation  procedure  and  examine  the  scalability, 
a  2D  and  3D  transonic  flows  are  calculated.  The  3rd  order  MUSCL  scheme  for  inviscid 
fluxes  and  2nd  order  central  differencing  scheme  for  the  viscous  terms  are  used  [136].  The 
governing  equations  are  solved  based  on  finite  volume  method. 

5.3.1  RAE2822  Transonic  Airfoil 

The  RAE2822  transonic  airfoil  is  calculated  to  examine  the  accuracy  of  the  CFD  solver 
and  parallel  computation  efficiency  for  2D  problems.  The  multi-block  grids  (Fig.  4.4)  are 
obtained  by  partitioning  a  single  block  O-grid  with  dimensions  of  257  x  56.  The  Reynolds 
number  is  6.5  x  106  based  on  the  chord  length.  The  Mach  number  is  0.729.  The  angle  of 
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attack  is  2.31°.  Fig.  5.8  shows  the  convergence  history  with  different  number  of  proces¬ 
sors  (blocks)  from  1  through  8.  The  convergence  rate  with  multiple  processors  is  slightly 
affected  by  the  approximate  implicit  treatment  at  the  sub-domain  boundaries.  Fig.  5.9 
presents  the  comparison  of  pressure  coefficients  between  the  experiment  and  computation 
with  1  and  8  processors.  The  computed  results  are  identical  and  agree  very  well  with  the 
experiment.  Excellent  speed  up  and  efficiency  for  parallel  computation  are  obtained  as 
shown  in  table  5.4  and  Fig.  5.10,  which  indicate  that  a  super-linear  scalability  is  achieved 
up  to  8  processors  with  the  sub-domain  grid  size  of  33  x  56. 


Figure  5.8:  Comparison  of  the  L2  residual  convergence  histories  of 
RAE2822  airfoil 
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Figure  5.9:  The  surface  distribution  of  pres¬ 
sure  coefficient  of  the  RAE2822  airfoil 


Figure  5.10:  Speedup  of  parallel  computation 
for  RAE2822  airfoil 
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Table  5.4:  The  parallel  computing  performance  for  2D  RAE2822  airfoil 


Number  of  nodes 

1 

2 

4 

8 

time  (sec/step) 

0.355 

0.161 

0.077 

0.039 

speed  up 

2.205 

4.610 

9.103 

Efficiency  (%) 

110.3 

115.3 

113.8 

Grid  size  (per  block) 

257  x  56 

129x56 

65  x56 

33x56 

5.3.2  Transonic  ONERA  M6  Wing 

This  case  is  to  examine  the  CFD  solver  accuracy  and  parallel  computation  efficiency  for  3D 
problems.  The  multi-block  grids  are  obtained  by  partitioning  a  single  block  O-H-grid  with 
the  dimensions  of  145  x  61  x  41(Fig.  4.11).  The  Mach  number  is  0.8395.  The  Reynolds 
number  is  1.97  x  107  based  on  the  averaged  chord.  The  angle  of  attack  is  0°. 

Fig.  5.11  presents  the  comparison  of  surface  pressure  distributions  between  the  exper¬ 
iment  and  computation  using  2-block  mesh  at  different  span-wise  sections.  The  location 
of  z/b  =  0.2  is  near  the  root  and  z/b  =  0.99  is  near  the  tip  of  the  wing.  The  computation 
results  agree  well  with  the  experimental  data.  Comparing  Fig.  5.11  and  Fig.  4.12  which 
uses  16-block  mesh,  the  computed  results  are  almost  identical. 

Fig.  5.12  shows  the  convergence  history  using  2  and  16  processors  (blocks).  The  con¬ 
vergence  rate  using  16  processors  is  also  slightly  affected  by  the  approximate  implicit  treat¬ 
ment  at  the  sub-domain  boundaries. 

Again,  excellent  speed  up  and  efficiency  for  the  parallel  computation  are  obtained  as 


shown  in  table  5.5  and  Fig.  5.13. 


Figure  5.11:  Surface  pressure  distribution  of  M6  wing  at  different  span- 
wise  locations 


Table  5.5:  The  parallel  computing  performance  for  3D  M6  wing 


Number  of  nodes 

1 

2 

4 

8 

16 

time  (sec/step) 

10.747 

4.912 

2.484 

1.298 

0.697 

speed  up 

2.188 

4.326 

8.280 

15.419 

Efficiency  (%) 

109.4 

108.2 

103.5 

96.37 

Grid  size  (per  block) 

145  x  61  x  41 

73x61  x  41 

37x61  x  41 

19x61  x  41 

10x61  x  41 
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Figure  5.12:  Comparison  of  the  maximum  residual  convergence  histories 
of  M6  wing 


Figure  5.13:  Speedup  of  the  parallel  computation  for  M6  wing 


Chapter  6 


Detached-Eddy  Simulation 


As  a  validation  of  the  DES  methodology,  a  circular  cylinder  flow  is  calculated.  The  5th 
order  WENO  finite  differencing  scheme  and  the  4th  order  central  differencing  scheme  de¬ 
scribed  in  Chapter3  with  the  LDE  Riemann  solver  and  implicit  time  marching  are  used. 

6.1  DES  of  a  Circular  Cylinder  Flow 

In  this  study,  the  flow  around  a  cylinder  at  a  Reynolds  number  of  3900  is  calculated  using 
DES.  The  Mach  number  is  0.2.  The  spanwise  length  is  nD,  where  D  is  the  cylinder  diam¬ 
eter.  The  dimensions  of  the  baseline  grid  are  (121  x  81  x  33)  (see  Fig.  6.1  and  Fig.  6.2). 
The  computation  is  conducted  on  an  MPI  based  computer  cluster  composed  of  200  Intel 
Xeon  5150  processors  with  the  floating  calculation  speed  of  2.66Ghz. 

A  non-dimensional  time  step  of  0.01  was  used  for  the  cases.  The  non-dimensional  time 
is  defined  as  t  —  Djv  .  The  computation  begins  with  a  uniform  flow  field.  All  the  results 
are  time-averaged  from  t  —  100  to  300. 

Fig.  6.3  shows  the  mean  pressure  coefficients  on  the  cylinder  surface.  For  Re  =  3900, 
only  the  coefficient  of  back  pressure  (Cp  at  0  =  180°)  is  available  .  The  computed  mean 
pressure  coefficient  agrees  very  well  with  the  experiment  at  0  <  6  <  60°.  The  present  result 
using  baseline  grid  is  better  than  the  FES  result  of  Kasliwal  et  al.  [143]  in  this  region.  In 
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Figure  6. 1 :  The  computational  grid  of  cylinder 


Figure  6.2:  Close-up  view  of  the  computational  grid 
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the  region  of  0  —  60°  ~  180°,  the  computed  pressure  lies  among  the  experimental  results. 


Figure  6.3:  Mean  pressure  coefficient  variation  on  the  surface  of  the  cylinder 

Fig.  6.4  is  the  averaged  mean  streamwise  velocity  on  the  centerline  in  the  wake  of  the 
cylinder.  The  present  result  agrees  better  with  the  experiment  [144]  than  those  of  LES 
[143, 145]  conducted  by  Kravchenko-Moin  and  Kasliwal  et  al. 

Fig.  6.5  through  Fig.  6.7  show  the  Reynolds  stress  components  located  at  x/D  =  1.54 
plane.  The  computed  streamwise  Reynolds  stress  ( uu  )  is  quite  symmetric  about  the  center 
line,  whereas  the  experiment  [146]  has  asymmetric  profile.  The  computed  ( u'u )  agrees 
well  with  the  experiment  except  it  does  not  reach  the  asymmetric  high  peak. 

The  computed  shear  Reynolds  stress  component  (mV)  in  Fig.  6.6  under-predicts  the 
amplitude  of  the  peaks  measured  in  the  experiment.  Fig.  6.7  also  shows  that  the  peak  of 
the  lateral  Reynolds  stress  (vV)  is  under  predicted.  However,  all  the  present  results  are 
significantly  better  than  the  LES  results  of  Rizzetta  et  al.  [147]  which  use  6-order  compact 
scheme  with  mesh  dimensions  of  199  x  197  x  53  as  shown  from  Fig.  6.5~  6.7. 
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Fig.  6.8  and  Fig.  6.9  plot  the  averaged  mean  streamwise  velocity  and  mean  crossflow 
velocity  at  three  streamwise  locations,  x/D  =  1.06,  x/D  —  1.54  and  x/D  =  2.02.  The 
present  results  agree  well  with  the  computed  results  of  Kravchenko-Moin  [145]  and  Kasli- 
wal  et  al.  [143]. 


y/D 


Figure  6.8:  Mean  streamwise  velocity  profiles  in  the  wake 


y/D 


Figure  6.9:  Mean  crossflow  velocity  profiles  in  the  wake 
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To  investigate  the  solution  sensitivity  to  spanwise  length,  the  cylinder  with  spanwise 
length  doubled  to  2 kD  is  calculated.  The  grid  density  is  the  same  as  the  baseline  grid. 
Fig.  6.10  through  Fig.  6.14  indicate  that  the  spanwise  length  has  only  a  small  effect  on 
the  computed  results.  The  computed  surface  pressure  is  nearly  identical  to  the  baseline 
results.  For  the  2 kD  spanwise  length,  the  computation  gives  slightly  lower  minimum  mean 
streamwise  velocity  in  the  wake  region.  There  is  also  little  difference  for  the  Reynolds 
stress  component  predicted  with  2 kD  spanwise  length. 

The  mesh  refinement  is  also  performed  in  this  study.  The  dimensions  of  the  baseline 
grid  are  increased  by  1.5  times  to  (181  x  121  x  49).  The  grid  is  divided  into  12  blocks 
for  parallel  computation.  The  mesh  refinement  has  a  significant  effect  on  the  computed 
results.  The  clear  difference  between  the  refined  and  baseline  mesh  is  that  the  refined  mesh 
has  larger  vortex  shedding  area.  Fig.  6.15  and  Fig.  6.16  show  the  contours  of  the  averaged 
mean  vorticity  magnitude  of  the  baseline  grid  and  refined  grid  respectively.  Both  display 
the  symmetry  of  the  mean  flow  field  after  a  long  time  average.  The  refined  grid  predicts  a 
larger  recirculation  zone  behind  the  cylinder. 

Fig.  6.10  shows  that  the  mean  pressure  distribution  is  raised  up  in  the  region  0  =  60°  ~ 
180°  and  matches  closer  with  the  experiment  of  Reynolds  number  3000.  The  streamwise 
velocity  distribution  is  shifted  away  from  the  measurement  value  as  shown  in  Fig.  6.11. 
The  computed  shear  Reynolds  stress  components  are  sharply  reduced  and  are  closer  to  the 
LES  results  of  Rizzetta  et  al. 

Fig.  6.17  shows  the  contours  of  the  instantaneous  vorticity  magnitude  at  I  =  300.  The 
refined  mesh  catches  more  small  scale  vortex  structures.  Fig.  6.18  shows  a  3D  instanta¬ 
neous  vorticity  magnitude  of  the  2nD  cylinder  at  t  =  300.  It  indicates  that  DES  resolves 
some  small  vortex  structures. 

In  summary,  for  the  baseline  grid  of  121  x81  x33  with  kD  spanwise  length,  the  com¬ 
puted  surface  pressure  and  velocity  in  the  wake  region  agree  well  with  the  experiment.  The 
computed  Reynolds  stress  are  also  in  good  agreement  with  the  experiment  except  that  the 
peak  values  are  some  what  under  predicted.  To  minimize  the  numerical  dissipation,  the  e 
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value  of  0.3  in  the  WENO  scheme  is  used. 

The  study  indicates  that  the  spanwise  length  of  kD  is  sufficient.  Doubling  the  spanwise 
length  yields  little  difference  of  the  results.  The  computation  of  mesh  refinement  indicates 
that  DES  is  significantly  affected  by  grid  size.  The  clear  difference  is  that  the  vortex  shed¬ 
ding  region  is  increased  when  the  mesh  is  refined.  The  increased  recirculation  zone  hence 
changes  the  mean  values  of  the  flow  field. 


Figure  6.10:  Mean  pressure  coefficient  variation  on  the  surface  of  the  cylinder 
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Figure  6.14:  Lateral  Reynolds  stress  at  x/D=1.54 
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Figure  6.15:  Contours  of  mean  vorticity  calculated  on  the  baseline  grid 
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Figure  6.16:  Contours  of  mean  vorticity  calculated  on  the  refined  grid 
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Figure  6.17:  Contours  of  instantaneous  vorticity  at  t=300T 
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Figure  6.18:  Overview  of  the  contours  of  instantaneous  vorticity  at  t=300T 


Chapter  7 


Validations  of  FSI  Model 


To  validate  the  high  order  fully  coupled  FSI  methodology  developed  in  this  research,  the 
vortex  induced  vibration  of  a  cylinder  and  a  forced  pitching  airfoil  are  simulated. 

7.1  Vortex-Induced  Oscillating  Cylinder 

In  this  section,  the  vortex- induced  oscillations  of  an  elastically  mounted  circular  cylinder  is 
computed  using  the  LDE  and  5th  order  WENO  scheme.  The  4th-order  fully  conservative 
central  differencing  is  employed  for  the  viscous  terms.  The  unsteady  laminar  Navier-Stokes 
equations  and  the  linear  structural  equation  are  fully  coupled  implicitly  via  successive  in¬ 
teraction  with  pseudo  time  stepping.  The  vortex-induced  oscillation  of  3D  cylinder  is  sim¬ 
ulated  for  the  first  time  in  this  research. 

7.1.1  2D  Simulation 

The  stationary  cylinder  is  simulated  first  to  provide  a  initial  flow  field  for  the  simulation  of 
vortex-induced  oscillating  cylinder. 
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7.1. 1.1  Stationary  Cylinder 

The  mesh  used  for  the  computation  of  stationary  cylinder  and  vortex-induced  oscillating 
cylinder  is  shown  in  Fig.  7.1.  The  dimensions  of  the  grid  are  121  x  81.  The  free- stream 
Mach  number  is  0.2.  The  Reynolds  number  based  on  the  diameter  of  the  cylinder  is  500. 
The  laminar  Navier-Stokes  equations  are  solved  because  of  the  low  Reynolds  number. 


X 

Figure  7.1:  The  computational  grid  of  cylinder 

The  time  history  of  the  computed  drag  and  lift  coefficients  is  shown  in  Fig.  7.2.  As 
shown  in  the  figure,  the  lift  oscillates  at  a  frequency  in  terms  of  the  Strouhal  number  Stq, 
which  is  in  a  good  agreement  with  the  experiments  of  Roshko  and  Goldstein  [148, 149] 
as  shown  in  table  7.1.  The  drag  coefficient  oscillates  with  twice  that  frequency,  Stq.  The 
computed  results  using  this  baseline  mesh  are  in  good  agreement  with  the  results  using 
refinement  mesh  in  Ref.  [33]  as  indicated  in  table  7.1.  The  reason  is  that  the  fifth-order 
WENO  and  newly  developed  LDE  schemes  have  better  accuracy.  Table  7.1  shows  the 
comparison. 

Fig.  7.3  and  Fig.  7.4  show  the  vorticity  contours  around  the  stationary  cylinder  with  the 
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Table  7.1:  Comparison  of  computed  results  with  the  experiments 


Mesh  dimension 

Stcd 

Stc, 

Cj  max 

Cdave 

120  x  80 

0.4496 

0.2248 

1.1822 

1.4658 

120x80  [33] 

0.4395 

0.2197 

1.181 

1.453 

200  x  120  [33] 

0.4516 

0.2246 

1.227 

1.484 

384  x  96  [28] 

0.4674 

0.2331 

1.149 

1.315 

(Roshko  1954  [148]) 

0.2075 

(Goldstein  1938  [149]) 

0.2066 

Figure  7.2:  Time  history  of  the  lift  and  drag  of  the  stationary  cylinder 

largest  and  smallest  lift  coefficients  respectively.  The  vortex  shedding  is  clearly  simulated 
in  these  figures. 

7. 1.1. 2  Vortex  Induced  Oscillating  Cylinder 

The  motion  of  the  elastically  mounted  cylinder  is  controlled  by  the  governing  equation 
described  in  section  2.3.1.  The  flow  conditions  for  this  oscillating  cylinder  are  the  same  as 
those  used  for  the  stationary  cylinder.  The  reduced  velocity  is  determined  by  St  number: 
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Figure  7.3:  Vorticity  contours  of  the  stationary  cylinder  with  the  largest  lift  coefficient 


Figure  7.4:  Vorticity  contours  of  the  stationary  cylinder  with  the  smallest  lift  coefficient 
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In  the  present  computation,  the  St  number  is  set  to  be  0.2,  corresponding  to  u  — 
1.5915.  The  mass  ratio,  fls  is  12.7324.  The  damping  ratio,  £  is  0.1583.  The  dimensionless 
physical  time  step  of  0.05  is  used. 

The  time  history  of  the  computed  lift  and  drag  coefficients  shown  in  Fig.  7.5  indicates 
that  the  averaged  drag  coefficient  is  larger  than  that  of  the  stationary  cylinder.  It  means  the 
motion  of  a  cylinder  enlarges  the  drag  coefficient.  The  amplitude  of  the  drag  coefficients 
is  also  enlarged  due  to  the  motion.  On  the  contrary,  the  amplitude  of  the  lift  coefficients  is 
decreased  because  of  the  motion. 

The  trajectory  of  the  central  point  of  the  cylinder  is  shown  in  Fig.  7.6.  The  trajectory  is 
similar  to  the  results  computed  by  [150]  and  [28]. 

Fig.  7.7  and  Fig.  7.8  show  the  vorticity  contours  around  the  oscillating  cylinder  with 
the  positive  and  negative  peak  lift  coefficients  respectively.  It  can  be  seen  that  the  vortexes 
keep  their  coherent  shedding  pattern  similar  to  the  one  simulated  around  the  stationary 
cylinder  because  the  displacement  is  small.  With  the  increase  of  the  displacement,  the 
vortex  shedding  pattern  will  become  more  and  more  irregular  [33]. 


Figure  7.5:  Time  history  of  the  lift  and  drag  of  the  oscillating  cylinder 
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Figure  7.6:  Time  history  of  the  displacement  of  the  oscillating  cylinder 
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Figure  7.7:  Vorticity  contours  of  the  oscillating  cylinder  with  the  largest  lift  coefficient 

7.1.2  3D  Simulation 

The  3D  cylinder  and  mesh  is  obtained  by  simply  extending  2D  cylinder  and  mesh  to 
one  diameter  long  in  spanwise  direction  (see  Fig.  7.9).  The  dimensions  of  the  grid  are 
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Figure  7.8:  Vorticity  contours  of  the  oscillating  cylinder  with  the  smallest  lift  coefficient 

121  x  81  x  11.  The  same  2D  structural  equations  given  in  section  2.3.1  is  employed  in 
3D  computation.  The  unsteady  lift  and  moment  coefficients  are  integrated  from  the  3d  un¬ 
steady  flow  fields.  To  compare  with  the  2D  results,  the  same  flow  and  structure  parameters 
are  use  in  the  3D  computation. 

7.1.2.1  Stationary  Cylinder 

Fig.7.10  plots  the  time  history  of  the  computed  lift  and  drag  coefficients.  The  compari¬ 
son  of  results  for  2D  and  3D  computation  is  shown  in  table  7.2.  They  indicate  that  the 
computed  Strouhal  number  of  3D  is  slightly  better  than  those  of  2D  when  compared  with 
the  experiments  [148, 149].  Both  lift  and  averaged  drag  coefficients  are  slightly  decreased 
compared  with  the  2D  results. 

The  vorticity  contours  around  the  stationary  cylinder  with  the  positive  and  negative 
peak  lift  coefficients  are  shown  in  Fig.  7.11  and  Fig.  7.12  respectively.  The  coherent  vortex 
shedding  pattern  is  the  same  as  the  2D  case. 
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Figure  7.9:  The  computational  grid  of  3D  cylinder 

7.1.2.2  Oscillating  Cylinder 

As  shown  in  Fig.  7.13,  the  amplitude  of  lift  coefficients  caused  by  the  oscillating  cylinder 
is  dramatically  decreased  compared  with  stationary  cylinder  results.  Different  from  the  2D 
case,  the  cylinder  motion  in  3D  has  a  smaller  effect  on  the  drag  coefficients.  Fig.  7.14  plots 
the  time  history  of  cylinder  displacement  and  the  amplitude  is  smaller  than  that  of  the  2D 
case. 

Fig.  7.15  and  Fig.  7.16  show  the  vorticity  contours  of  the  oscillating  cylinder  with  the 
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Table  7.2:  Comparison  of  results  for  2D  and  3D  computation 


Stcd 

Stc, 

Cj  max 

Cdave 

2D  computation 

0.4496 

0.2248 

1.1822 

1.4658 

3D  computation 

0.4422 

0.2211 

1.0585 

1.4294 

(Roshko  1954  [148]) 

0.2075 

(Goldstein  1938  [149]) 

0.2066 

Figure  7.10:  Time  history  of  the  lift  and  drag  of  the  stationary  cylinder 


high  and  low  peak  lift  coefficients  respectively.  Same  as  the  2D  case,  the  vortexes  keep 
their  coherent  shedding  pattern  similar  to  the  one  obtained  for  the  stationary  cylinder. 
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Figure  7.11:  Vorticity  contours  of  the  stationary  cylinder  with  the  largest  lift  coefficient 


Figure  7.12:  Vorticity  contours  of  the  stationary  cylinder  with  the  smallest  lift  coefficient 
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Figure  7.13:  Time  history  of  the  lift  and  drag  of  the  oscillating  cylinder 


Figure  7.14:  Time  history  of  the  displacement  of  the  oscillating  cylinder 
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Figure  7.15:  Vorticity  contours  of  the  oscillating  cylinder  with  the  largest  lift  coefficient 


Figure  7.16:  Vorticity  contours  of  the  oscillating  cylinder  with  the  smallest  lift  coefficient 


114 


7.2  Forced  Pitching  Vibration  of  NACA  64A010  Airfoil 

As  a  validation  case  of  DES  for  fluid- structural  interaction,  a  forced  pitching  airfoil,  NACA64A010 
is  calculated.  The  NACA64A010  airfoil  is  selected  as  the  validation  case  because  the  ex¬ 
perimental  data  are  available.  The  forced  pitching  airfoil  is  simulated  first  using  2D  RANS. 

7.2.1  2D  Simulation  Using  RANS 

For  this  transonic  airfoil,  an  O-type  grid  is  generated  with  the  dimensions  of  281  x  66  (see 
Fig.  7.17). 


Figure  7.17:  The  computational  grid  of  NACA64A010  airfoil 


The  airfoil  is  forced  in  pitch  about  its  quarter  chord  sinusoidally.  The  airfoil  oscillation 
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is  defined  by  a  function  of  the  time  dependent  variation  of  its  AoA, 

a(t)  =  ocq  + aAsm(cot)  (7.1) 

where  a  ( t )  is  the  time  dependent  AoA.  ccq  is  the  mean  of  the  oscillating  angle,  CCa  is  the 
amplitude  of  the  oscillating  angle,  co  is  the  angular  frequency,  which  is  directly  related  to 
the  reduced  frequency 


where  C  is  the  chord  of  the  airfoil,  and  £/«>  is  the  free-stream  velocity. 

To  be  consistent  with  the  experiment,  the  following  primary  parameters  are  employed 
in  the  unsteady  calculation:  (Xq  =  0 ,  0Ca  =  1-01°,  Reynolds  number  (based  on  chord),  Re  — 
1.256  x  107,  free-stream  Mach  number,  M =  0.8,  reduced  frequency,  kc  —  0.202.  The 
computation  begins  with  a  uniform  flow  field  of  free  stream  at  0°  AoA.  The  dimensionless 
time  step  is  At  =  0.05. 

Fig.  7.18  shows  the  lift  coefficients  varying  with  the  AoA  after  the  flow  field  reaches 
its  temporally  periodic  solution.  The  computed  lift  coefficients  agree  well  with  the  experi¬ 
ment  [151].  Fig.  7.19  shows  the  moment  coefficients  varying  with  the  AoA.  The  agreement 
of  the  moment  coefficient  is  not  as  good  as  that  of  the  lift  coefficient.  However,  the  agree¬ 
ment  in  the  current  results  is  better  than  the  recent  result  computed  by  McMullen  et  al. 
in  2002  [152].  The  discrepancy  between  the  computation  and  the  experiment  in  the  mo¬ 
ment  coefficient  may  be  caused  by  the  inadequacy  of  the  shock/turbulence  boundary  layer 
interaction,  which  may  not  predict  the  surface  friction  accurately. 


Figure  7.18:  Comparison  of  computed  lift  coefficient  with  experimental  data  for  the  forced 
pitching  airfoil 


Figure  7.19:  Comparison  of  computed  moment  coefficient  with  experimental  data  for  the 
forced  pitching  airfoil 


Chapter  8 


DES  of  Fluid-Structural  Interaction 


This  chapter  simulates  NLR7301  airfoil  limit  cycle  oscillation  (LCO)  caused  by  fluid- 
structural  interaction  (FSI)  using  DES  method.  As  a  validation  case,  the  forced  pitching 
NACA64A010  airfoil  is  simulated  firstly.  Then  the  2D  NLR7301  airfoil  limit  cycle  oscil¬ 
lation  (LCO)  is  simulated  using  RANS  method.  Finally,  the  3D  NLR7301  airfoil  LCO  is 
simulated  using  DES  method.  The  low  diffusion  E-CUSP  (LDE)  scheme  with  5th  order 
weighted  essentially  non-oscillatory  scheme  (WENO)  is  employed  to  calculate  the  inviscid 
fluxes.  The  fully  conservative  4th  order  central  differencing  is  used  for  the  viscous  terms. 
The  fully  coupled  fluid-structural  interaction  model  is  employed  in  all  the  cases. 

8.1  Validation  for  Forced  Pitching  Airfoil 

A  DES  validation  is  conducted  to  simulate  the  fluid- structural  interaction  of  the  forced 
pitching  NACA64A010  airfoil.  Since  DES  must  be  3D  for  LES  of  turbulence,  the  airfoil 
is  extended  in  spanwise  direction.  An  O-H-type  is  generated  with  the  dimensions  of  281  x 
66  x  49  (see  Fig.  8.1).  The  spanwise  length  extended  is  4  times  of  the  chord  length  of  the 
airfoil.  The  forced  pitching  equation  and  parameters  are  the  same  as  those  used  in  the  2D 
computation  described  in  section  7.2.1. 

Fig.  8.2  shows  the  lift  coefficients  varying  with  the  AoA  after  the  3D  flow  field  reaches 
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Figure  8.1:  The  computational  grid  of  NACA64A010  airfoil 


its  temporally  periodic  solution.  The  computed  lift  coefficients  using  DES  agree  very  well 
with  the  experiment. 

Fig.  8.3  shows  the  moment  coefficients  varying  with  the  AoA.  The  computed  moment 
coefficients  using  DES  have  a  deviation  compared  with  the  experiment.  The  discrepancy 
between  the  computation  and  the  experiment  in  the  moment  coefficient  may  be  caused  by 
the  inadequacy  of  the  shock/turbulence  boundary  layer  interaction,  which  may  not  predict 
the  surface  friction  accurately 

Since  there  are  no  separations  in  the  flow  field,  the  lift  and  moment  using  DES  can 
be  considered  as  the  spanwise-averaged  values.  Thus,  the  lift  and  moment  using  DES  are 
almost  the  same  as  those  using  RANS  in  section  7.2.1. 

8.2  Limit  Cycle  Oscillations  of  NLR7301  Airfoil  Using  RANS 

In  this  section,  the  fluid-structural  interaction  of  the  elastically  mounted  NLR7301  airfoil 
is  simulated  by  using  RANS  with  S-A  one  equation  model.  The  RANS  simulation  is  to 
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Figure  8.2:  Comparison  of  computed  lift  coefficient  with  experimental  data  for  the  forced 
pitching  airfoil 

search  the  initial  condition  to  match  the  experimental  LCO.  The  initial  condition  is  then 
used  for  DES  to  save  CPU  time. 

The  case  simulated  is  the  test  case  No. 77  [40,41, 153]  of  NLR7301  airfoil.  The  chord 
length  of  the  airfoil  is  0.3 m  and  the  mean  angle  of  attack  is  1.28°.  The  experimental  con¬ 
ditions  are  at  a  free- stream  Mach  number  of  0.768  and  a  Reynolds  number  of  1.727  x  106 
based  on  the  chord  length.  The  experiment  was  conducted  at  a  total  pressure  of  pt  =  0.45 
bar  and  a  dynamic  pressure  of  pdyn  =  0. 126  bar. 

The  elastically  mounted  NLR7301  airfoil  has  two  degree  of  freedom,  plunge  and  pitch. 
The  structural  motion  equation  are  given  in  section  2.3.2. 

The  non-dimensional  structural  parameters  used  for  the  computation  of  fluid-structural 
interaction  are  summarized  in  table  8.1 

For  numerical  simulation,  there  are  several  initial  conditions  affecting  LCO:  free  stream 
Mach  number,  initial  AoA  and  «o.  The  present  research  investigates  the  effects  of  these 
initial  conditions  and  finds  the  optimum  values  of  them  that  generate  the  LCO  matching 
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Figure  8.3:  Comparison  of  computed  moment  coefficient  with  experimental  data  for  the 
forced  pitching  airfoil 


Table  8.1:  The  non-dimensional  structural  parameters  used  for  computation 


Mach  number 

r  a 

<5/, 

(Oh 

0.768 

0.0484 

0.197 

0.0041 

0.0073 

942 

0.31988 

0.24306 

0.753 

0.0484 

0.197 

0.0041 

0.0073 

942 

0.32625 

0.24790 

best  with  the  experiment.  Two  different  procedures  are  used  to  simulate  the  LCO. 

Procedure  1 

The  first  procedure  follows  the  criterion  used  by  Weber  et  al.  [39]  and  Tang  et  al.  [2],  in 
which,  both  free  stream  Mach  number  and  initial  AoA  are  adjusted  to  match  the  computed 
steady  state  surface  pressure  distribution  with  experiment  as  much  as  possible.  The  ob¬ 
tained  Mach  number  and  AoA  are  then  used  in  the  LCO  simulation  by  adjusting  the  «o 
until  the  computed  time-averaged  AoA  is  close  to  the  obtained  AoA  of  the  steady  state 
computation. 


121 


Procedure  2 

The  second  procedure  is  developed  in  the  present  research,  in  which  the  free  stream  Mach 
number  is  fixed  and  the  oco  as  well  as  initial  AoA  both  are  iterated  to  match  the  experimental 
LCO  amplitudes.  The  resulted  time-averaged  lift  and  moment  are  also  taken  into  account 
to  compare  with  the  experiment. 

All  simulations  are  conducted  on  an  MPI  based  computer  cluster  with  parallel  compu¬ 
tation.  The  parallel  computation  is  performed  by  the  high  efficiency  parallel  computing 
algorithm  described  in  Chapter  5  [134]. 

8.2.1  Steady  State  Flow  Computation 

The  baseline  grid  is  a  single  block  O-type  grid  with  grid  dimensions  of  193  x  97  and  is 
equally  partitioned  into  16  sub-blocks  with  8  blocks  in  the  circumferential  direction  and 
two  blocks  in  radial  direction  as  shown  in  Fig.  8.4.  The  dimensions  of  each  sub-domain  is 
25  x  49. 


Figure  8.4:  The  computational  grid  of  NLR7301  airfoil 
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The  steady  state  computations  are  conducted  first  to  search  for  the  flow  conditions  that 
match  the  computed  pressure  distribution  best  with  the  experiment.  Both  the  widely  used 
Mach  number  of  0.753  for  CFD  [2, 39]  and  the  experimental  Mach  number  of  0.768  are 
simulated.  It  is  found  that  the  steady  state  surface  pressure  agrees  best  with  experiment 
at  AoA  =  —0.45°  for  Mach  0.753  and  at  AoA  —  —0.2°  for  Mach  0.768.  However,  both 
cases  have  some  discrepancy  with  the  experiment.  Fig.  8.5  shows  the  surface  pressure 
distribution  for  the  two  cases.  For  the  case  of  Mach  number  0.753,  the  shock  location  at 
suction  surface  agrees  better  with  the  experiment,  whereas  for  the  Mach  0.768,  the  shock 
location  at  pressure  surface  agrees  better.  Overall,  the  computed  case  with  the  experimental 
Mach  number  of  0.768  is  closer  to  the  experiment  upstream  and  downstream  of  the  shocks. 

The  mesh  refinement  is  performed  for  the  steady  state  case  at  Mach  0.768  to  confirm 
that  the  baseline  mesh  is  sufficient  to  be  used  for  unsteady  LCO  simulation.  The  baseline 
mesh  is  refined  in  both  directions  with  the  mesh  size  increased  by  4  times  to  385  x  193. 
As  shown  in  Fig.  8.5,  the  computed  surface  pressure  distributions  between  the  baseline  and 
refined  mesh  have  little  difference,  except  that  the  refined  mesh  has  sharper  shock  profile 
due  to  the  denser  mesh. 

8.2.2  2D  LCO  Simulation 

The  first  series  of  LCO  search  use  the  procedure  1  at  Mach  =  0.753  and  AoA  =  —0.45°, 
which  are  the  conditions  used  by  the  research  groups  of  Weber  et  al.  and  Tang  et  al.  [2,39]. 
The  LCO  computation  is  conducted  by  adjusting  «o.  the  off-wind  value  of  a  to  make  the 
time-averaged  AoA  agree  with  AoA  —  —0.45°.  The  LCO  computation  is  conducted  using 
two  different  initial  flow  fields  to  investigate  the  effects  of  initial  flow  and  perturbation. 
One  initial  field  is  the  solution  of  the  steady  state  computation.  The  other  initial  field  is  the 
uniform  free  stream  flow.  Both  LCO  computations  start  at  AoA  =  0°.  The  dimensionless 
physical  time  step  of  0.01  is  used  which  is  defined  as  tc  =  t/(c/uoo).  t  is  the  physical 
time.  Fig.  8.6  and  8.7  show  the  computed  LCO  amplitudes  at  Mach  number  of  0.753 
and  (Xq  of  0.25°.  Both  initial  fields  predict  the  final  LCOs  with  about  the  same  amplitudes 
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Figure  8.5:  Pressure  coefficient  distribution  on  the  surface  of  NLR7301  airfoil 

even  though  the  transition  period  is  different.  Fig.  8.8  and  8.9  plot  the  lift  and  moment 
coefficients  respectively.  The  convergence  history  within  a  typical  physical  time  step  is 
shown  in  Fig.  8.10. 

The  second  series  of  LCO  search  use  the  procedure  2  at  the  Mach  0.753.  The  LCO  is 
conducted  by  iterating  «o  and  initial  AoA  to  match  the  measured  LCO  amplitudes  as  close 
as  possible. 

Table  8.2  lists  all  the  trail  iteration  cases  at  Mach  0.753.  Case  A  shows  the  computed 
results  using  procedure  1.  Case  B  to  Case  F  show  the  computed  results  with  the  initial  field 
setup  as  uniform  free  stream  at  the  different  AoA  and  Cfo  using  procedure  2. 

It  can  be  seen  from  Table  8.2  that  the  Case  A  matches  the  lift  and  moment  coefficients 
best  with  the  experiment  among  all  the  cases.  However,  the  predicted  LCO  amplitudes  are 
an  order  of  magnitude  higher  than  the  experiment,  just  like  all  other  cases. 

Table  8.2  indicates  that  both  AoA  and  Cfo  have  influence  on  the  amplitudes  of  LCO. 
However,  the  computation  is  not  able  to  match  the  amplitudes  by  adjusting  AoA  and  ocq 
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Figure  8.6:  Pitch  motion  predicted  by  RANS  (M  =  0.753, AoA  =  —0.45°) 


Figure  8.7:  Plunge  motion  predicted  by  RANS  ( M  =  0.753,  AoA  =  —0.45°) 
due  to  the  Mach  number  that  is  different  from  the  experiment. 


Fig.  8.11  shows  the  contours  of  Mach  number  labeled  from  (a)-(j)  with  the  interval  of 
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Figure  8.8:  Lift  coefficient  predicted  by  RANS  (M  —  0.753, AoA  =  —0.45°) 


Figure  8.9:  Moment  coefficient  predicted  by  RANS  ( M  —  0.753, AoA  =  —0.45°) 


1/10  of  a  cycle  for  Case  E,  which  has  the  smallest  amplitudes  in  all  the  cases  at  the  Mach 
number  0.753.  Fig.  8.12  and  Fig.  8.13  plot  the  corresponding  positions  of  (a)  to  (j)  in 
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Table  8.2:  Comparison  of  the  test  cases  at  M  =  0.753 


NLR7301  LCO 

Initial  AoA 

«o 

Lift  Coef. 

Moment  Coef. 

h(mm) 

a° 

Case  A 

0.0 

0.25 

0.2318 

-0.0800 

10.134 

3.1942 

Case  B 

0.0 

0.60 

0.2944 

-0.0744 

9.6768 

2.9796 

CaseC 

0.0 

0.75 

0.3180 

-0.0756 

8.6943 

2.6921 

Case  D 

0.0 

0.85 

0.3365 

-0.0770 

7.8386 

2.4352 

Case  E 

0.0 

0.95 

0.3548 

-0.0790 

6.9359 

2.1566 

Case  L 

0.05 

0.85 

0.3455 

-0.0784 

7.3517 

2.2863 

Experiment 

0.272 

-0.082 

0.75 

0.20 

a  cycle  for  pitching  and  plunging  movement  respectively.  The  phase  difference  between 
pitching  and  plunging  movement  is  168°.  The  experimental  data  is  176°.  Under  this  large 
amplitude  LCO,  the  shock  location  and  strength  vary  significantly.  At  position  (a),  the  AoA 
is  maximum.  There  is  only  one  strong  shock  on  suction  surface.  With  the  AoA  decreased 
from  (b)  to  (e),  the  suction  surface  shock  is  weakened.  At  position  (e),  the  AoA  is  the 
minimum.  The  double  shock  pattern  is  formed  on  suction  surface  and  the  boundary  layer  on 
suction  surface  is  the  thinest.  This  is  because  the  shock  strength  is  the  weakest  with  double 
shock.  While  the  suction  surface  shock  is  weakened,  a  shock  appears  on  pressure  surface. 
After  position  (e),  the  AoA  is  increased.  The  shock  on  suction  surface  is  strengthened 
with  a  single  shock.  The  boundary  layer  on  suction  surface  becomes  thick  again  due  to 
the  strong  shock/turbulent  boundary  layer  interaction.  From  position  (h)  to  (]),  while  the 
suction  surface  shock  becomes  stronger,  the  pressure  surface  shock  disappears.  As  shown 
in  Fig.  8.14,  when  the  AoA  is  maximum,  there  is  a  small  boundary  layer  separation  on  the 
trailing  edge. 

The  third  series  of  LCO  search  still  use  the  procedure  2,  but  at  the  experimental  Mach 
number  of  0.768.  The  dimensionless  physical  time  step  is  the  same  as  that  used  at  Mach 
0.753. 

Table  8.3  lists  the  iteration  cases  at  Mach  0.768.  The  computed  averaged  lift,  moment 
coefficients,  frequencies  and  amplitudes  of  Case  E  all  agree  excellently  with  the  experi¬ 
ment.  Fig.  8.15,  8.16,  8.17  and  8.18  show  the  computed  LCO  amplitudes  of  pitch,  plunge, 
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lift  and  moment  coefficient  at  Mach  number  of  0.768  for  Case  E.  Compared  with  results 
predicted  at  Mach  =  0.753  as  shown  in  Fig.  8.6,  8.7,  8.8  and  8.9,  it  can  be  seen  that  the 
predicted  amplitudes  at  Mach  —  0.768  are  more  than  one  order  of  magnitude  smaller.  This 
means  that  the  Mach  number  has  a  significant  effect  on  the  amplitudes  of  plunge  and  pitch¬ 
ing  oscillation.  The  reason  may  be  that  the  different  Mach  number  causes  different  shock 
strength,  different  shock/boundary  layer  interaction  patterns  and  hence  different  unsteady 
non-linear  forcing  and  moment. 
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Figure  8.15:  Pitch  motion  predicted  by  RANS  (M  =  0.768) 


Totally  685,000  physical  time  steps  are  performed  in  the  computation  of  case  E  with 
tc  =  6850.  The  stabilized  ECO  period  in  the  simulation  is  up  to  2250/ r  and  88  cycles.  The 
convergence  history  within  a  typical  physical  time  step  is  shown  in  Fig.  8.19.  It  indicates 
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Figure  8.16:  Plunge  motion  predicted  by  RANS  ( M  =  0.768) 


that  only  about  7  pseudo  time  steps  are  needed  to  reach  the  convergence  criteria  of  10-6  in 
this  computation.  Fig.  8.20  shows  the  contours  of  Mach  number  at  one  cycle  for  case  E. 
Fig.  8.21  and  Fig.  8.22  plot  the  corresponding  positions  of  (a)  to  (j)  in  a  cycle  for  pitching 
and  plunge  movement  respectively.  The  computed  phase  difference  between  pitching  and 
plunge  movement  is  172°,  which  agrees  very  well  with  the  experimental  phase  difference 
of  176°.  Different  from  the  LCO  computation  at  Mach  0.753,  there  is  no  separation  in  the 
simulated  flow  field  at  Mach  0.768  as  shown  in  Fig.  8.23.  The  reason  is  that  the  very  small 
amplitudes  of  the  LCO  captured  at  Mach  0.768  does  not  cause  large  AoA  variation. 

Table  8.4  summarizes  the  computed  LCO  amplitudes  and  frequencies  at  different  con¬ 
ditions  compared  with  the  experimental  results  [40].  At  Mach  —  0.753,  the  present  com- 
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Figure  8.17:  Lift  coefficient  predicted  by  RANS  ( M  =  0.768) 

puted  results  are  comparable  to  those  of  Weber  et  al.  and  Tang  et  al.  However,  at  Mach 
0.768  of  the  experimental  condition,  both  the  predicted  plunge  and  pitching  amplitudes 
agree  excellently  with  the  experiment,  whereas  the  previous  results  predicted  by  other  re¬ 
searchers  [2,39]  at  Mach  —  0.753  are  more  than  one  order  of  magnitude  higher.  This  is  the 
first  time  that  a  numerical  simulation  of  NLR7301  airfoil  LCO  matches  the  experiment. 

In  general,  the  prediction  accuracy  of  Case  E  at  Mach  0.768  is  on  the  same  order  of 
the  experiment  measurement  uncertainty.  The  only  primary  difference  from  the  experi¬ 
ment  for  Case  E  is  that  the  «o  used  in  the  simulation  is  0.75°  whereas  the  experimental 
value  is  1.28°.  The  «o  only  affects  the  initial  moment  imposed  on  the  elastic  system  and 
remains  as  a  constant  in  the  whole  LCO  process.  Such  difference  may  be  attributed  to  the 
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time 

Figure  8.18:  Moment  coefficient  predicted  by  RANS  ( M  =  0.768) 

uncertainty  of  the  experiment  and  numerical  simulation,  and  the  sensitive  nature  of  LCO  to 
initial  perturbations,  which  are  difficult  if  not  impossible  to  be  made  the  same  between  the 
experiment  and  numerical  simulation. 

Note  that  the  final  LCO  plunge  amplitude  is  about  2.7/1000  of  the  chord  length  and  the 
pitching  amplitude  is  0.2287°.  These  are  some  very  small  values.  The  accurate  resolution 
of  such  small  scale  vibration  motion  without  being  damped  out  in  the  long  time  calculation 
may  be  attributed  to  the  high  order  low  diffusion  numerical  schemes  and  the  fully  coupled 
FSI  model  employed  in  this  research. 

It  needs  to  point  out  an  important  phenomenon  that  the  LCO  amplitudes  are  dependent 
on  the  initial  flow  fields.  The  results  in  Table  8.3  are  computed  using  the  initial  field  set 
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Table  8.3:  Comparison  of  the  test  cases  at  M  =  0.768 


NLR7301  LCO 

Ini.  AoA 

«o 

Ci 

Cm 

h(mm ) 

a° 

f(Hz) 

Case  A 

0.0 

0.68 

0.2610 

-0.0796 

1.6509 

0.4632 

33.35 

Case  B 

0.0 

0.70 

0.2758 

-0.0807 

4.2453 

1.2349 

33.36 

CaseC 

0.0 

0.75 

0.2729 

-0.0805 

1.2617 

0.3524 

33.38 

Case  D 

-0.033 

0.75 

0.2673 

-0.0799 

1.4351 

0.4015 

33.39 

Case  E 

0.05 

0.75 

0.2803 

-0.0816 

0.8192 

0.2287 

33.36 

Experiment 

1.28 

0.272 

-0.082 

0.75 

0.20 

32.74 

Table  8.4:  LCO  comparison  of  computation  and  experiment 


NLR7301  LCO 

Mach 

h(mm ) 

err  (mm) 

a° 

err 

Wz) 

err  (Hz) 

Present 

0.768 

0.8192 

0.069 

0.2287° 

0.0287° 

33.36 

0.62 

Present 

0.753 

10.134 

9.384 

3.1942° 

2.9942° 

33.49 

0.75 

Weber  (2001)  [39] 

0.753 

10.5 

9.75 

4.09° 

3.89° 

33.42 

0.68 

Tang  (2003)  [2] 

0.753 

8.99 

8.24 

3.17° 

2.97° 

34.3 

1.56 

Experiment  [40] 

0.768 

0.75 

0.20° 

32.74 

equal  to  the  uniform  free-stream  everywhere.  If  a  converged  steady  state  solution  is  used 
as  the  initial  field,  the  LCO  amplitudes  may  be  very  different  with  significantly  greater 
magnitude.  The  different  «o  and  initial  AoA  also  setup  the  initial  lift  and  moment  to  cer¬ 
tain  values.  In  other  words,  different  initial  perturbation  may  generate  very  different  LCO 
solutions.  This  appears  to  be  the  bifurcation  phenomenon  due  to  the  non-linear  aerody¬ 
namic  loading  of  lift  and  moment,  which  are  caused  by  the  pattern  of  shock  wave/turbulent 
boundary  layer  interaction.  If  we  can  understand  the  systematic  relationship  between  the 
LCO  amplitude  and  initial  perturbation,  it  may  be  possible  to  control  the  LCO  and  mitigate 
or  prevent  it.  Even  though  the  LCO  amplitudes  are  very  different  as  shown  in  Table  8.3  and 
8.4,  the  computed  frequency  varies  little  and  agrees  well  with  the  experimental  value. 
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8.3  Limit  Cycle  Oscillations  of  NLR7301  Airfoil  Using  DES 

Based  on  the  experience  of  the  2D  LCO  simulation  using  RANS,  the  3D  DES  of  LCO  for 
the  same  experimental  conditions  is  conducted.  The  same  structural  equations  (2.44)-(2.46) 
are  used  with  fully  coupled  fluid-structural  interaction  procedure.  The  unsteady  lift  and 
moment  coefficients  are  integrated  from  the  3d  unsteady  flow  fields.  The  computation  grid 
is  composed  of  24  blocks  partitioned  from  a  single  block  O-H-grid  with  the  dimension  of 
193  x  97  x  33  (see  Fig.  8.24).  This  simulation  is  very  CPU  intensive  and  will  be  impossible 
without  parallel  computation.  The  span  wise  length  is  3.33  times  of  the  airfoil  chord,  which 
is  the  same  as  that  used  in  the  experiment.  The  dimensionless  physical  time  interval  is  0.01 . 
The  Procedure  2  is  adopted  in  the  3D  LCO  simulation.  Based  on  the  experience  of  2D  LCO 
simulation,  the  computation  is  conducted  with  the  experimental  Mach  number  0.768  only. 
Three  cases  are  performed  with  different  initial  AoA  and  do. 

Case  A:  AoA  =  0.05°,  ao  =  0.75° 

In  this  case,  AoA  =  0.05°,  ao  =  0.75°  are  used  in  the  computation  which  are  the  same  as 
the  2D  Case  E  at  Mach  =  0.768.  The  predicted  amplitudes  of  plunge  and  pitch  modes  are 
shown  in  Fig.  8.25  and  Fig.  8.26.  Different  from  the  corresponding  2D  case,  the  ampli¬ 
tudes  of  the  3D  case  are  4  times  higher  than  the  experiment.  Fig.  8.27  and  Fig.  8.28  plot 
the  variation  of  lift  and  moment  coefficients  with  the  time  after  the  LCO  is  formed.  As 
indicted  in  table  8.5,  the  predicted  time-averaged  lift  coefficient  is  greater  than  those  of 
the  corresponding  2D  case.  The  predicted  time-averaged  moment  coefficient  is  lower  than 
those  of  the  2D  case. 


Table  8.5:  Comparison  of  LCO  using  RANS  and  DES 


NLR7301  LCO 

Ini.  AoA 

«o 

h  (mm) 

a° 

c. 

Cm 

/(Hz) 

DES 

0.05 

0.75 

2.8424 

0.8066 

0.2934 

-0.0831 

33.41 

RANS 

0.05 

0.75 

0.8192 

0.2287 

0.2803 

-0.0816 

33.36 

Experiment 

1.28 

0.75 

0.20 

0.272 

-0.082 

32.74 

141 


Figure  8.24:  The  3D  computational  grid  of  NLR7301  airfoil 
Case  B:  AoA  =  0°,  «o  =  0.75° 

In  this  case,  the  initial  AoA  is  decreased  to  zero  to  observe  the  variation  of  the  LCO.  The 
predicted  amplitudes  of  plunge  and  pitch  modes  are  shown  in  Fig.  8.29  and  Fig.  8.30. 
Fig.  8.31  and  Fig.  8.32  plot  the  variation  of  lift  and  moment  coefficients  with  the  time  after 
the  LCO  is  formed.  As  shown  in  table  8.6,  the  amplitudes  of  the  LCO  are  greater  than  those 
of  case  A.  Compared  with  case  A,  the  predicted  time-averaged  lift  coefficient  is  decreased, 
but  the  predicted  time-averaged  moment  coefficient  is  increased.  The  absolute  value  of  the 
moment  coefficient  is  decreased. 
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Figure  8.27:  Lift  coefficient  predicted  by  DES  ( AoA  —  0.05°,  ao  =  0.75°) 


Figure  8.28:  Moment  coefficient  predicted  by  DES  ( AoA  =  0.05°,  «o  =  0.75°) 
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Figure  8.29:  Pitch  motion  predicted  by  DES  ( AoA  =  0.0°,  Oq  =  0.75°) 


Figure  8.30:  Plunge  motion  predicted  by  DES  ( AoA  =  0.0°,  (Xq  =  0.75°) 
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Case  C:  AoA  =  0°,  aQ  =  0.85° 

In  this  case,  the  ao  is  increased  to  0.85°,  which  is  closer  to  the  experimental  data  1 .28°.  The 
initial  AoA  is  kept  the  same  as  the  case  B.  The  predicted  amplitudes  of  plunge  and  pitch 
modes  are  shown  in  Fig.  8.33  and  Fig.  8.34.  Fig.  8.35  and  Fig.  8.36  plot  the  variation  of  lift 
and  moment  coefficients  with  the  time  after  the  LCO  is  formed.  As  shown  in  table  8.6,  the 
amplitudes  of  the  LCO  are  significantly  decreased  compared  with  case  A  and  case  B.  They 
are  much  closer  to  the  experiment  than  case  A  and  case  B.  The  predicted  time-averaged  lift 
coefficient  is  larger  than  case  A.  The  predicted  time-averaged  moment  coefficient  is  lower 
than  case  A. 


Figure  8.33:  Pitch  motion  predicted  by  DES  ( AoA  =  0.0°,  Oo  =  0.85°) 

The  instantaneous  vorticity  magnitudes  at  three  spanwise  locations,  which  are  located 
at  the  left  most,  right  most  and  middle  sections  respectively,  are  shown  in  Fig.  8.37  when 
the  time  dependent  AoA  reaches  its  maximum.  The  vorticity  contours  show  a  strong  shock 
on  suction  surface  and  a  weak  shock  on  pressure  surface  which  are  similar  to  those  demon¬ 
strated  in  Fig.  8.20.  Some  small  vortex  structures  due  to  the  shock  wave/boundary  layer 
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Figure  8.36:  Moment  coefficient  predicted  by  DES  ( AoA  —  0.0°,  Oo  =  0.85°) 
interaction  are  resolved  by  the  DES. 

Case  D:  AoA  =  -0.015°,  «o  =  0.88° 

In  this  case,  the  initial  AoA  is  decreased  to  —0.015°.  The  Oo  is  increased  to  0.88°.  The 
predicted  amplitudes  of  plunge  and  pitch  modes  are  shown  in  Fig.  8.38  and  Fig.  8.39. 
Fig.  8.40  and  Fig.  8.41  plot  the  variation  of  lift  and  moment  coefficients  with  the  time  after 
the  ECO  is  formed.  As  shown  in  table  8.6,  the  amplitudes  of  the  ECO  agree  excellently 
with  the  experiment.  The  ocq  used  in  the  simulation  is  closer  to  the  experiment  than  the  best 
2D  case.  The  predicted  time-averaged  lift  and  moment  coefficients  have  slight  differences 
from  the  experiment.  Same  as  the  2D  cases,  the  computed  frequency  varies  little  and  agrees 
well  with  the  experimental  value.  In  summary,  the  computed  results  of  this  3D  ECO  case 
match  the  experiment  very  well,  better  than  the  2D  case. 
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Figure  8.37:  Instantaneous  vorticity  magnitude  predicted  by  DES  ( AoA  =  0.0°,  ao  =  0.85°) 


Figure  8.40:  Lift  coeffl 
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Figure  8.41:  Moment  coefficient  predicted  1 
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Table  8.6:  Comparison  of  LCO  for  different  initial  AoA  and  ccq  at  M  =  0.768  using  DES 


NLR7301  LCO 

Ini.  AoA 

«o° 

h  (mm) 

a° 

Q 

Cm 

f(Hz) 

Case  A 

0.05 

0.75 

2.8424 

0.8066 

0.2934 

-0.0831 

33.41 

Case  B 

0.0 

0.75 

3.4313 

0.9871 

0.2780 

-0.0788 

33.39 

CaseC 

0.0 

0.85 

1.0776 

0.2920 

0.3044 

-0.0850 

33.51 

Case  D 

-0.015 

0.88 

0.7344 

0.1980 

0.3073 

-0.0860 

33.53 

Experiment 

1.28 

0.75 

0.20 

0.272 

-0.082 

32.74 

Chapter  9 


Conclusions 


This  research  developed  an  efficient  and  accurate  methodology  to  resolve  flow  non-linearity 
of  fluid- structural  interaction.  A  numerical  strategy  to  apply  the  detached-eddy  simulation 
with  a  fully  coupled  fluid-structural  interaction  model  is  established  for  the  first  time.  The 
following  novel  numerical  algorithms  are  created:  a  general  sub-domain  boundary  mapping 
procedure  for  parallel  computation  to  reduce  wall  clock  simulation  time,  an  efficient  and 
low  diffusion  E-CUSP  (LDE)  scheme  used  as  a  Riemann  solver  to  resolve  discontinuities 
with  minimal  numerical  dissipation,  and  an  implicit  high  order  accuracy  WENO  scheme  to 
capture  shock  waves. 

The  Detached-Eddy  Simulation  is  based  on  the  model  proposed  by  Spalart  in  1997. 
Near  solid  walls  within  wall  boundary  layers,  the  Reynolds  averaged  Navier-Stokes  (RANS) 
equations  are  solved.  Outside  of  the  wall  boundary  layers,  the  3D  filtered  compressible 
Navier-Stokes  equations  are  solved  based  on  large  eddy  simulation(LES).  The  Spalart- 
Allmaras  one  equation  turbulence  model  is  solved  to  provide  the  Reynolds  stresses  in  the 
RANS  region  and  provide  the  subgrid  scale  stresses  in  the  LES  region.  The  validation  of 
the  DES  strategy  is  performed  by  simulating  the  flow  field  of  a  cylinder.  The  computed 
results  agree  well  with  the  experiment  except  that  the  peak  values  of  Reynolds  stress  are 
some  what  under  predicted.  The  computation  of  mesh  refinement  indicates  that  DES  is 
significantly  affected  by  grid  size  due  to  the  modeled  stress  depletion  (MSD)  effect.  The 
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implementation  of  DDES  as  the  next  step  is  expected  to  remove  this  problem. 

The  numerical  experiments  show  that  the  developed  LDE  scheme  is  able  to  capture 
crisp  shock  profiles  and  exact  contact  surfaces  with  the  improved  5th  order  finite  differenc¬ 
ing  WENO  scheme.  The  fully  conservative  4th  order  finite  central  differencing  schemes 
used  for  the  viscous  terms  are  verified  to  be  efficient  and  accurate.  The  LDE  scheme  is  more 
efficient  than  the  Roe  scheme  when  coupled  with  the  S-A  one  equation  model.  The  extra 
equation  of  the  S-A  turbulence  model  changes  the  Jacobian  of  the  Roe  scheme,  weakens 
the  diagonal  dominance,  reduces  the  maximum  CFL  number  permitted  by  the  Roe  scheme, 
and  hence  decreases  the  convergence  rate.  Both  the  LDE  and  the  Roe  scheme  predict  good 
results  when  compared  with  the  experiment  in  the  validation  cases.  The  computed  results 
of  the  LDE  scheme  agree  very  well  with  those  of  Roe  scheme.  The  numerical  experiments 
also  show  that  the  strategy  of  using  conservative  finite  differencing  discretization  to  an  ex¬ 
isting  finite  volume  code  is  achieved  without  changing  the  structure  of  the  code.  Hence, 
the  code  runs  with  high  efficiency  when  using  high  order  WENO  scheme.  For  the  time 
accurate  unsteady  simulation,  the  Jameson’s  dual  time  step  method  with  the  unfactored 
Gauss-Seidel  relaxation  iteration  shows  very  good  performance  in  this  research. 

A  general  sub-domain  boundary  mapping  procedure  is  developed  for  arbitrary  topol¬ 
ogy  multi -block  structured  grids  with  grid  points  matched  on  sub-domain  boundaries.  The 
interface  of  two  adjacent  blocks  is  uniquely  defined  according  to  each  local  mesh  index 
system  (MIS)  which  is  specified  independently.  A  pack/unpack  procedure  based  on  the 
definition  of  the  interface  is  developed  to  exchange  the  data  in  a  ID  array  to  minimize  data 
communication.  A  secure  send/receive  procedure  is  employed  to  remove  the  possibility  of 
blocked  communication  and  achieve  optimum  parallel  computation  efficiency.  Two  terms, 
“Order’  and  “Orientation” ,  are  introduced  as  the  logics  defining  the  relationship  of  adja¬ 
cent  blocks.  The  domain  partitioning  treatment  of  the  implicit  matrices  is  to  simply  discard 
the  comer  matrices  so  that  the  implicit  Gauss-Seidel  iteration  can  be  implemented  within 
each  subdomain.  The  numerical  experiments  indicate  that  the  effect  of  this  implicit  treat¬ 
ment  on  the  convergence  rate  is  small.  The  message  passing  interface  (MPI)  protocol  is 
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used  for  the  data  communication.  The  code  is  portable  to  any  platform  as  long  as  the  MPI 
is  available.  The  numerical  experiments  including  2D  and  3D  RANS  and  DES  simulations 
show  that  the  general  mapping  procedure  developed  in  this  research  is  robust  and  have  high 
scalability 

The  DES  solver  with  fully  coupled  fluid- structural  interaction  methodology  is  validated 
with  vortex  induced  vibration  of  a  cylinder  and  a  transonic  forced  pitching  airfoil.  For  the 
cylinder,  the  laminar  Navier-Stokes  equations  are  solved  due  to  the  low  Reynolds  number. 
The  3D  effects  are  observed  in  both  stationary  and  oscillating  cylinder  simulation  because 
of  the  flow  separations  behind  the  cylinder.  For  the  transonic  forced  pitching  airfoil  DES 
computation,  there  is  no  flow  separation  in  the  flow  field.  The  DES  results  agree  well  with 
the  RANS  results.  These  two  cases  indicate  that  the  DES  has  more  effects  on  the  flow  fields 
with  separation. 

The  code  is  then  used  to  simulate  the  limited  cycle  oscillation  of  NLR7301  airfoil  in  2D 
RANS  mode  with  the  S-A  one  equation  model  and  3D  DES  mode.  At  experimental  Mach 
number  of  0.768,  the  computed  LCO  amplitudes  and  frequency  are  in  excellent  agreement 
with  experiment  by  adjusting  AoA  and  CX().  The  time  averaged  lift  and  moment  coefficients 
also  match  the  experiment  very  well  using  RANS  model.  The  free-stream  Mach  number 
has  the  major  effect  on  the  amplitudes  of  LCO  due  to  different  shock/boundary  layer  in¬ 
teraction  patterns.  The  initial  flow  field  or  initial  perturbation  has  a  strong  influence  on 
the  amplitudes  of  LCO.  The  prediction  accuracy  is  on  the  same  order  of  the  experiment 
measurement  uncertainty.  The  only  primary  difference  from  the  experiment  is  that  the  ao 
used  in  the  simulation  is  0.75°  for  RANS  model  and  0.88°  for  DES  model,  whereas  the 
experimental  value  is  1.28°.  The  (Xq  only  affects  the  initial  moment  imposed  on  the  elas¬ 
tic  system  and  remains  as  a  constant  in  the  whole  LCO  process.  Such  difference  may  be 
attributed  to  the  uncertainty  of  the  experiment  and  numerical  simulation,  and  the  sensitive 
nature  of  LCO  to  initial  perturbations,  which  are  difficult  if  not  impossible  to  be  made  the 
same  between  the  experiment  and  numerical  simulation. 

This  research  appears  to  be  the  first  time  that  a  numerical  simulation  of  NLR7301  airfoil 
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LCO  matches  the  experiment.  It  may  be  attributed  to  the  high  order  low  diffusion  numerical 
schemes,  the  fully  coupled  FSI  model,  and  the  turbulence  model  used  in  this  research. 

The  simulations  of  LCO  also  confirm  some  of  the  experimental  observations  and  an¬ 
swers  some  important  questions.  First,  The  LCOs  with  the  small  relative  amplitude  is 
captured  with  unbounded  flows  in  the  numerical  simulation.  This  means  they  should  not 
be  the  artifacts  of  the  wind-tunnel  experiment  and  most  likely  are  the  factual  phenomenon. 
Second,  the  co-existence  of  multiple  LCOs  at  constant  flow  conditions  is  confirmed  in  our 
simulation.  The  reason  that  other  numerical  simulations  only  capture  the  LCOs  with  large 
amplitudes  may  be  due  to  their  high  numerical  dissipation  that  either  smears  out  the  small 
amplitude  LCO  or  is  only  able  to  resolve  the  large  amplitudes  LCOs.  Third,  the  numerical 
simulation  of  this  research  confirms  that  the  wall  boundary  layer  transition  from  laminar 
to  turbulent  does  not  have  a  large  effect  on  LCOs  at  high  Reynolds  number  because  our 
simulation  assumes  that  the  boundary  layer  is  fully  turbulent  from  the  airfoil  leading  edge. 
Fourth,  the  simulation  confirms  that  the  wind  tunnel  wall  interference  with  or  without  per¬ 
forated  test  section  does  not  have  much  effect  on  LCOs  because  the  present  simulation  uses 
the  unbounded  flow  condition  with  no  wind  tunnel  wall  effect  at  all.  Fifth,  the  numerically 
captured  LCO  at  Mach  0.768  is  not  accompanied  with  any  flow  separation  due  to  the  very 
small  amplitude.  This  may  modify  the  hypothesis  that  the  LCOs  are  caused  by  the  non¬ 
linearity  of  flow  separation  induced  by  shock/boundary  layer  interaction.  In  other  words, 
the  nonlinearity  of  shock/boundary  layer  interaction  with  no  flow  separation  is  sufficient 
to  trigger  a  LCO.  This  may  also  make  reduced  numerical  models  feasible  to  capture  LCOs 
with  reasonable  accuracy. 

In  conclusion,  the  numerical  strategy  of  the  high  order  DES  with  fully  coupled  FSI 
model  and  parallel  computing  developed  in  this  research  is  demonstrated  to  have  high  ac¬ 
curacy,  robustness,  and  efficiency. 


Chapter  10 


Future  Work 


Although  a  good  progress  has  been  achieved  in  the  present  research,  the  following  are  the 
areas  that  need  to  be  further  improved. 


Multi-grid  Method 

Computation  efficiency  is  always  an  important  issue  in  the  numerical  simulation.  High 
order  CFD  schemes  are  generally  more  CPU  time  consuming  because  of  the  large  amount 
of  mesh  points  used  and  the  complicated  numerical  algorithms.  The  multi-grid  method 
has  been  shown  to  be  a  effective  convergence  acceleration  method.  It  is  expected  that 
the  multi-grid  method  with  the  implicit  Gauss-Seidel  relaxation  may  further  improve  the 
computation  convergence  and  reduce  the  CPU  time. 


Detached-Eddy  Simulation 

Detached-eddy  simulation  is  a  relatively  new  method.  Numerical  experiments  show  that 
it  is  affected  by  the  grid  due  to  the  MSD.  To  overcome  the  MSD  problem  and  make  the 
DES  limiter  independent  of  grid  spacing,  Spalart  suggested  a  modification  to  the  original 
DES97  model  in  2006  [63],  referred  to  as  Delayed  DES(DDES).  The  DDES  model  has 
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demonstrated  excellent  agreement  with  experiment  and  a  significant  improvement  over  the 
DES97  for  the  tested  cases.  The  implementation  of  DDES  is  expected  to  remove  the  MSD 
problem. 


Wall  Boundary  Condition 

At  the  current  work,  it  is  observed  that  sometimes  the  high  order  wall  boundary  condition 
(BC)  are  not  as  robust  as  the  2nd  order  wall  BC.  Further  work  is  needed  to  make  the  high 
order  wall  BC  more  robust. 


Control  of  Limit  Cycle  Oscillation 

By  using  different  initial  flow  condition,  various  factors  affecting  LCO  have  been  investi¬ 
gated.  A  group  of  parameters  have  been  found  sensitive  to  the  LCO  of  NLR7301  airfoil. 
This  make  it  possible  to  control  the  formation  of  LCO  by  adjusting  the  initial  flow  field, 
which  could  be  a  new  interesting  research  area. 


Appendix  A 


Nonlinear  equations  of  airfoil  vibration 
motion 


A  solid  flexibly  supported  airfoil  is  shown  in  Fig.  A.l,  A. 2  and  A. 3.  In  Fig.  A.l,  the 
position  of  the  elastic  axis  (EO),  the  center  of  gravity  (T)  and  the  airfoil  chord  (c)  are 
sketched.  The  airfoil  can  be  vertically  displaced  and  rotated  about  EO.  Fig.  A. 2  shows  the 
elastic  support  of  the  airfoil  on  translational  and  rotational  springs. 


Figure  A.l:  The  position  of  the  elastic  axis  (EO),  the  centre  of  gravity  (T)  and  the  chord  (c) 

The  pressure  and  viscous  forces  acting  on  the  vibrating  airfoil  are  determined  by  the 
components  of  the  stress  tensor  and  result  in  the  lift  force  L(t)  and  the  torsional  moment 
M(t).  The  airfoil  in  neutral  and  deformed  positions  is  shown  in  Fig.  A. 3.  Based  on  the 
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neutral  position,  the  horizontal  and  vertical  displacements  of  any  point  on  the  airfoil  chord 
can  be  expressed  as 


u  =  x(  1  — cosa)  ,w  —  h  +  xsina 


(A.l) 


respectively.  Here  x  denotes  the  local  coordinates  measured  along  the  airfoil  chord  c  from 
the  elastic  axis.  The  kinetic  energy  Ek  of  the  airfoil  has  the  form 


Ek  =JA 


i  ( du\ 
dr  J  ^  \  dr  J 


Ps  (jc)  dx 


(A.2) 


where  ps  denotes  the 
Eq.  (A.2)  we  obtain 


=u 


( h+xacosa )  +( xasina)  ’  ps(x)dx 


density  of  the  airfoil  per  unit  length.  With  further  rearrangement  of 


Ek  =  \h2  fc Ps  (x)  dx  +  hacosa  fcxps(x)dx+ 

\a2cos2oc  fc.x2ps  (x)  dx  +  \a2sin2a  fcx2ps  (x)  dx  (A.3) 

=  \h2m  +  hacosaSa  +  \dc2la 

where 

m  —  Jcps  (x)  dx  is  the  mass  of  the  airfoil, 

Sa  —  fcxps  (x)  dx  is  the  static  moment  around  the  elastic  axis  EO, 

I a  —  Jcx2ps  (x)  dx  is  the  inertia  moment  around  the  elastic  axis  EO, 

The  potential  energy  V  of  the  airfoil  is 


V=\khh2+l-kaa2  (A.4) 

where  ki,  and  ka  are  the  bending  stiffness  and  torsional  stiffness  respectively. 

Kinetic  and  potential  energy  have  to  satisfy  the  Lagrange  equations 

'p-**  +  pL  =  Qj  (A.5) 

dt  dqj  dqj  dqj 


where  qj  (j=l,2)  are  generalized  coordinates,  i.e.  h  and  a  in  this  case,  and  Qj  are  general 
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ized  forces,  i.e.  the  aerodynamic  force  L(t)  and  the  moment  M(t). 


Figure  A. 2:  The  elastic  support  of  the  airfoil  on  translational  and  rotational  springs 
Thus,  for  j=l,2  we  have 


(A.6) 


[ hm  +  acosaSa]  +  kj,h  =  — L(t ) 

2f  [hcosaSa  +  «/«]  +  hasinocSa  +kh  =  M(t) 

Differentiation  with  respect  to  time  in  Eq.  (A.6)  yields  the  nonlinear  equations  of  mo 
tion  of  the  airfoil 


mh+Saacosa  —  Sacc2sina  +khh  =  —L(t) 

Sahcosoc  +  Iacc  +  kaa  =  M(t ) 

For  small  values  of  the  angle  a  and  of  its  derivative  a  (i.e.  sin  a  ~  a,  cos  a  ~ 
0),  Eq.  (A. 7)  yield  the  well  known  linearized  system 


(A.7) 


1,  da  « 


mh  +  Sa0C  +  kfrh  =  —L(t) 
Sah  T7a(i  -\-kaC(,  =  M{t) 


(A.8) 


Including  viscous  damping  terms  leads  to  the  governing  equations 


Appendix  B 


Derivation  of  Spalart-Allmaras 
Equation  (Eq.(2.23)) 


For  the  left  hand  side  of  Eq.(2.22), 


Pw  =P{w  +  v'V* 

=  pt*+pV-Vv 
=  ^v_v^  +  v.(pi>F 


dv 


_  dpv 
~  dt 


=  ^  +  v. 


pvV 


+  V-  (^pvV^j  -  v 


-  vV  •  ypV 

t  +  v-(pv 


For  the  right  hand  side  of  Eq.(2.22), 


£[V-((v  +  v)Vv)+cw(Vv)2] 

=  £V((v  +  v)Vv)  +  gcw(Vv)2 
=  V  •  [§(v  +  v) Vv]  -  i(v  +  v)Vv  ■  Vp  +  £c,2(Vv)2 

Let  Sv  be  the  source  term: 

Sv  =  pcbiSv  (1  —  fti)  —  p  (cwifw-^kfa)  {%)2~ 
i(v  +  v)Vv  ■  Vp  +  £q,2(Vv)2  +pfn  (A q)2 


(B.l) 


(B.2) 


(B.3) 
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The  conservative  form  of  Eq.(2.22)  can  be  written  as 


dpv 

~di 


+ V-  pvV  = V 


—  (v  +  v)Vv 

a 


+  S\ 


(B.4) 


The  dimensionless  flow  variables  in  the  governing  equations  are  defined  as  the  follow¬ 


ing, 


X  = 

X 

L’ 

y  = 

y 

L’ 

2  = 

Z 

L’ 

ll 

v 

w 

u  — 

uZ' 

V  — 

EoT’ 

w  = 

z  Z7oT’ 

n  — 

p 

= 

V 

t  — 

t 

P  — 

Poo’ 

8 

8 

L/Uoo 

where  the  bar  denotes  dimensionless  variable,  the  free  stream  conditions  are  denoted 
by  °o,  L  is  the  reference  length  used  in  the  Reynolds  number  Re/,, 


ReL 


Poot/ooL 

f-icc 


(B.5) 


For  the  left  hand  side  of  Eq.(B.4), 


i§i  +  v.(pw 


—  dPl  i 

~  dt  + 


dpuV  .  dpvv  .  dpwV 


_  Poo  Ucx 


_  PooUoo 


dx 


+ 


dpv  , 
dt  + 

dpv 


dy 

dpiiv 

dx 

/ 


+ 


+ 


dz 
dpvV 
dy 
\ 


+ 


dpwV 

dz 


T  +  V-  PVV 


(B.6) 


where, 


-  -zd  -?d 
V~‘d7x  +  JTy 


For  the  right  hand  side  of  Eq.(B.4), 


V- 


—  (v  +  v)Vv 
La 


-(v  +  T)V9 

a 


(B.7) 


The  first  term  of  Sv: 


165 


pcbiSv  (1  —  fti)  =  pcb\  (1  —  fa)  (^^+L2p  Jjifo)  *  (B-8) 

The  second  term  of  Sv: 


The  third  term  of  Sv: 


^(v  +  v)Vv- Vp 


i 

L2Poo  o 


(v  +  v)  Vv-Vp 


The  fourth  term  of  Sv: 


(B.10) 


The  fifth  term  of  Sv: 


~Q,2(Vv)2 


p 

T2Poo  CJ 


cb2  (Vv) 


2 


(B  .11) 


P/m  (Ap)2  -  Poof/2p/i  (Ag)2  (B.12) 

The  Eq.(B.4)  then  can  be  written  as  the  following, 


M  +  v- 

dt  + 


pvE 


=  J_v* 

Re  V 


^(v  +  v)Vv  +pcb\  (1  ~  ftl)  (5+  )  v— 


ifep  (c-i/w  -  %fa)  (I)  -  £ i  ( v  +  *)  V*  •  Vp+ 

i§cw(W)2+^p/fl(A# 


(B.13) 


For  convenience,  we  drop  the  bar  used  in  Eq.(B.13).  Following  the  deriving  process  of 
Hu’s  PhD  thesis  [136],  Eq.(B.13)  can  be  transformed  from  Cartesian  coordinate  system  to 
the  generalized  coordinate  system  as  the  following 


djpV  ,  dpvU  ,  dpvV  ,  dpvW  _ 
dt  +  dt,  +  dp  +  dt;  ~ 

1  ( 5§(v+v)(l»Vv)  ,  3§(v+v)(m.Vv)  ,  d§(v+v)(n«Vv)  ,  i  c  ^ 
Re  {  d %  +  dt 7  +  dt;  +  J^v ) 


(B.14) 
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It  is  the  Eq.(2.23). 
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